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ABSTRACT 


An analysis was made to determine the natural frequencies and mode 
shapes of ring- and/or stringer-stiffened noncircular cylinders with 
arbitrary end conditions. The method of analysis used and the results 
of the analysis are presented in Volume I of this report (Reference 1). 
Volume II contains the computer program and the user instructions for the 
program. Sample input and output is presented in the appendices. 
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A THEORETICAL ANALYSIS OF THE FREE VIBRATIONS OF 
RING- AND/OR STRINGER- STIFFENED ELLIPTICAL 
CYLINDERS WITH ARBITRARY END CONDITIONS , 

VOLUME II - USERS MANUAL FOR 
COMPUTER PROGRAM 

By Donald E. Boyd, C. K. P. Rao and 
Robert L. Brugh 

/ 

INTRODUCTION 

The free vibration characteristics of ring- and/or stringer- 
stiffened circular and noncircular cylindrical shells are of interest to 
designers of flight and marine structures. Frequently, fuselages of 
flight structures and hulls of submarines have noncircular cross-section £ 
due either to special internal storage requirements or to imperfections 
occurring during manufacture. The method of analysis developed in 
Volume I of this report (Reference 1) is capable of evaluating the free- 
vibrational characteristics of ring- and/or stringer-stiffened "singly" 
symmetric noncircular cylinders with arbitrary end conditions. In this 
analysis, the stiffeners are treated as discrete elements. The 
stiffeners may be arbitrarily located and all stiffeners need not possess 
the same geometric and material properties; however, the stiffeners are 
assumed to be uniform along their axes. The analysis considers the 
extension and flexure of the shell and extension, torsion, and flexure 
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about both cross-section axes of the stiffeners. The stringers may have 
non symmetric cross-sections but the rings are assumed to have "singly" 
symmetric cross-sections. The rotary inertia of the shell is neglected. 

Based on this method of analysis, a computer program was developed. 
Using this program, a comparative study was made using known solutions 
for circular and noncircular, unstiffened and stiffened cylinders with 
various end conditions. Results of this study are presented in 
Reference 1. The limitations of the program and instructions fo^ using 
the program are discussed in the following paragraphs of this report. 



V PROGRAM LIMITATIONS 

A flow chart of the main program is given in Appendix A. The main 
program and some of its subroutines, as listed in Appendix B of this 
report, are written in single precision for use on the CDC Model 6600 
Computer. The mass and stiffness matrices for the entire shell struc- 
ture. are generated in the program. This program uses the subroutine 
EIGENP (2) to determine the eigenvalues and eigenvectors of the problem. 
A dictionary of the variables used in the main program is presented in 
Appendix C. 

The computer program has the following limitations : 

A. Shell 

1. Constant thickness 

2. Isotropic material properties 

B. Stringers 

1. Maximum number; 16 

2. Maximum kinds; 1 

3. Uniform along length 

C . Rings 

1. Maximum number; 11 

2. Maximum kinds; 2 

3. Uniform around circumference 

4. "Singly" symmetric about z-axis 


3 


i 



D. Number of terms 

The maximum, number of terms in general must satisfy the 
following equation. 

3 (n terms used) (m terms used) £ 90 

For a specific case, refer to the equations given in the 
computer program for determining the value of MN3. MN3 must 
be less than or equal to 90. 

The limitations on the program may be made less restrictive by 
increasing the appropriate dimensions in the dimension statement of 
the main program. 



USER INSTRUCTIONS 

In addition to the main program arid subroutines listed, the user 
must supply three function subroutines. The first, FUNCTION RSHL(T), 
defines the radius of curvature [r] of the shell as a function of the 
0 coordinate. FUNCTION RRRT(T) defines the first derivative with respect 
to 0 of the reciprocal of the radius f I t) | . . The third, FUNCTION 

•• X >e : . 

RSHLT(T) , defines the first derivative with respect to 0 of the radius 
[R, Q ] . As an example, page 79(a) presents the subroutines written for 
an elliptical cylinder having a specific major (A) and minor (B) axis. 

The input data for the program is prepared according to Appendix D. 
The input data is divided into the following four categories: (1) 

general data; (2) shell data; (3) stringer data; and (4) ring data. The 
general and shell data are required for all computer runs. The program 
in its current state will solve problems with the following boundary 
conditions: free-free, clamped-free, freely supported, and clamped- 

clamped. The input variables are defined at the beginning of the program 
listing. 

The other two categories (stringer and ring data) are needed only 
when the shell structure is stiffened by rings and/or stringers. A set 
of stringer and/or ring data will be required for each kind of ring 
and/or stringer used to stiffen the shell. 

A computer output for an example problem is presented in Appendix E. 
The example problem has both stringer and ring stiffening. It should be 
noted that all input data is given on the printout. The first page gives 
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the general Information and shell data. The second and third pages give 
the stringer and ring data, respectively. The stringer and ring data 
pages will appear in the printout only when the stiffening is used in the 
problem. Other printout options may be selected such that the stiffness 
matrix, the mass matrix and the eigenvectors may be printed out. 

An input listing for this same problem is presented in Appendix F. 
Cards one through four give the general information. The shell data is 
on card five. Cards six through thirteen and fourteen through nineteen 
give the stringer and ring data, respectively. The twentieth card is 
the first card of the general information of the second problem. The 
integer "one" (1) punched in column 80 of thiB card Indicates it is the 
last card of the data set. It should be noted that there is no limit to 
the number of problems which can be solved in each run. 
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APPENDIX A 

FLOW CHART OF THE MAIN PROGRAM 
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♦ ♦ 
♦ REFERENCE ♦ 

*■ 
* 


•FREE VIBRATIONAL ANALYSIS OF STIFFENED OR UNST I FF EN ED 
CIRCULAR OR NONCIRCULAR CYLINDERS WITH ARBITRARY END 

CONDITIONS* 


LANGUAGE USED 
DIGITAL MACHINE 
PROGRAMMER 


FORTRAN IV 
IBM 360/65 

C. K. PANDURANGA RAO 
GRADUATE RESEARCH ASSOCIATE 
SCHOOL OF MECHANICAL AND 
AEROSPACE ENGINEERING 
OKLAHOMA STATE UNIVERSITY 
STILLWATER, OKLAHOMA 7407* 

AUGUST 30, 1971 + 

♦ 
♦ 
♦ 

♦ THIS PROGRAM COMPUTES BOTH THE SYMMETRIC AND ANTISYMMETRIC ♦ 

♦ FREQUENCIES AND THE CORRESPONDING EIGENVECTORS OF STIFFENED OR ♦ 

♦ UNSTIFFENED CIRCULAR OR NONCIRCULAR CYLINDERS WITH ARBITRARY END ♦ 

♦ CONDITIONS. THE RAYLEIGH-RIT l METHOD IS USED TO GENERATE THE ♦ 

♦ STIFFNESS AND MASS MATRICES. ♦ 

♦ ♦ 


DATE OF COMPLETION 


DESCRIPTION OF THE PROGRAM 
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DESCRIPTION OF THE PARAMETERS 
INPUT PARAMETERS 


BCR 

NQUIT 

NG 


KG 

LL 

NL 

KK 

NK 

MM IN 

MMAX 

MSA 


NMIN 

NMAX 

NSA 


NEO 


IR 


KKR 


NAME OF THE BOUNDARY CONDITION 

1 IN THE 80 TH COLUMN OF A BLANK CARD AT THE END 
OF THE DATA SETS TO SIGNIFY THE END OF DATA SETS 
OROER OF THE GAUSSIAN QUADRATURE. NG HAS TO BE 
ANY ONE OF THE FOLLOWING NUMBERS 3,4,5,6,7,8,9,10, 
16, AND 32. 

NUMBER OF CIRCUMFERENTIAL INTERVALS INTO WHICH THE 

LIMITS OF INTEGRATION ARE DIVIDED 

TOTAL NUMBER OF STRINGERS 

NUMBER OF KINDS OF STRINGERS 

TOTAL NUMBER OF RINGS 

NUMBER OF KINDS OF RINGS 

STARTING VALUE OF MIN THE ASSUMED DISPL SERIES 
FINAL VALUE OF M IN THE ASSUMED DISPL SERIES 

0 WHEN ONLY EVEN M VALUES ARE CONSIDERED 

1 WHEN ONLY ODD M VALUES ARE CONSIDERED 

2 WHEN BOTH EVEN ANO ODD VALUES OF M ARE CONSIDERED 

STARTING VALUE OF N IN THE ASSUMED DISPL SERIES 
FINAL VALUE OF N IN THE ASSUMED OISPL SERIES 

0 WHEN COMPUTING THE SYMMETRIC MODE SHAPES WITH 
RESPECT TO THE VERTICAL AXIS OF THE CROSS-SECTION 

1 WHEN COMPUTING THE ANTISYMMETRIC MODE SHAPES WITH 
RESPECT TO THE VERTICAL AXIS OF THE CROSS-SECTION 

0 WHEN ONLY EVEN N VALUES ARE CONSIDERED 

1 WHEN ONLY ODD N VALUES ARE CONSIDERED 

2 WHEN BOTH EVEN AND ODD VALUES OF N ARE CONSIDERED 

0 WHEN THE CROSS-SECTION OF THE SHELL IS CIRCULAR 

1 WHEN THE CROSS-SECTION OF THE SHELL IS 
NONCIRCULAR 

NUMBER OF THE KINDS OF RINGS WHICH HAVE DIFFERENT 


♦ 

♦ 

♦ 

♦ 

♦ 

♦ 
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♦ 

*■ 
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♦ 
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♦ 
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*■ 
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♦ 
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4 KWK 
♦ 

4 

4 

4 NWM 

4 

4 ' 

4 

4 NWEV 

4 

4 

4 

4 TITLE1 
4 TITLE2 
4 PC 
4 EC 
4 XNU 
4 H 
4 AA 

4 NNL (L I 
4 

4 TIL, II 

4 

4 psm 

4 ES( L ) 

4 ASIL) 

4 Z 1S( L ) 

4 

4 Z2SIL) 

4 

4 Y1S(L) 

4 

4 

4 Y2S(LI 
♦ 

4 Z IS ( L ) 

4 

4 YISIL) 

■f 

4 

4 YZIS(L) 
♦ 

4 

4 GJSIU 
♦ 

4 NNKIKI 
4 

4 RX ( K, I) 
4 

4 PR(K) 

4 ERIK I 
4 AR(K) 

4 E1RIK) 

♦ 

4 E2RIKJ 

4 

4 ZIRIKI 

4 

4 

4 XIRIK) 


CENTROIDS 4 

0 WHEN THE STIFFNESS MATRIX IS NOT TO BE PRINTED 4 

1 WHEN THE STIFFNESS MATRIX IS TO BE PRINTED 4 

2 WHEN THE STIFFNESS MATRIX IS TO BE PRINTEO ANC 4 

PUNCHED OUT ON THE CARDS 4 

0 WHEN THE MASS MATRIX IS NOT TO BE PRINTED 4 

1 WHEN THE NASS MATRIX I S TO BE PRINTED 4 

2 WHEN THE MASS MATRIX IS TO BE PRINTED AND 4 

PUNCHED OUT ON THE CARDS 4 

0 WHEN THE EIGENVECTOR MATRIX IS NOT TO BE PRINTED 4 

1 WHEN THE EIGENVECTOR MATRIX IS TO BE PRINTED 4 

2 WHEN THE EIGENVECTOR MATRIX IS TO BE PRINTED AND 4 

PUNCHED OUT ON THE CARDS 4 

THE TITLE OF THE PROBLEM 4 

THE TITLE OF THE PROBLEM ICONTINUEDI 4 

THE MASS DENSITY OF THE SHELL 4 

THE YOUNG'S MODULUS OF THE SHELL ♦ 

THE POISSON'S RATIO OF THE SHELL 4 

THE THICKNESS OF THE SHELL 4 

LONGITUDINAL LENGTH OF THE SHELL 4 

NUMBER OF STRINGERS WHICH HAVE THE L TH SET OF 4 

PROPERTIES 4 

LIST OF THETA VALUES UN DEGREES) AT WHICH THE L TH 4 
SET OF STRINGERS ARE LOCATEO 4 

THE MASS DENSITY OF THE L TH SET OF STRINGERS 4 

THE YOUNG'S MOOULUS OF THE L TH SET OF STRINGERS 4 

CROSS-SECTIONAL AREA OF THE L TH SET OF STRINGERS ♦ 
THE Z-OI STANCE OF THE SHEAR CENTER OF THE L TH SET 4 
OF STRINGERS FROM THE SHELL'S MIDDLE SURFACE 4 

THE Z-OISTANCE OF THE CENTROIO OF THE L TH SET OF 4 
STRINGERS FROM THEIR SHEAR CENTER 4 

THE Y-OISTANCE OF THE SHEAR CENTER OF THE L TH SET 4 
OF STRINGERS FROM THE Z-AXIS PASSING THROUGH THEIR 4 
POINTS OF ATTACHMENT 4 

THE Y-DIST ANCE OF THE CENTROID OF THE L TH SET OF 4 
STRINGERS FROM THEIR SHEAR CENTER 4 

THE MOMENT OF INERTIA OF THE CROSS-SECTION OF THE 4 
L TH SET OF STRINGERS ABOUT THE Z-AXIS PASSING 4 

THROUGH THEIR CENTROID 4 

THE MOMENT OF INERTIA OF THE CROSS-SECTION OF THE 4 
L TH SET OF STRINGERS ABOUT THE Y-AXIS PASSING 4 

THROUGH THEIR CENTROID 4 

THE PROOUCT INERTIA OF THE CROSS-SECTION OF THE 4 

L TH SET OF STRINGERS ABOUT Y- AND Z-AXES PASSING 4 
THROUGH THEIR CENTROID 4 

THE TORSIONAL STIFFNESS OF THE L TH SET OF 4 

STRINGERS 4 

NUMBER OF RINGS WHICH HAVE THE K TH SET OF RING 4 

PROPERTIES 4 

LIST OF X-POSITIONS OF THE RINGS WITH K TH SET OF 4 
RING PROPERTIES 4 

THE MASS OENSITV OF THE K TH SET OF RINGS 4 

THE YOUNG'S MODULUS OF THE K TH SET OF RINGS 4 

THE CROSS-SECTIONAL AREA OF THE K TH SET OF RINGS 4 

THE Z-DI STANCE OF THE SHEAR CENTER OF THE K TH SET 4 

OF RINGS FROM THE MIDDLE SURFACE OF THE SHELL 4 

THE Z-OISTANCE OF THE CENTROID OF THE K TH SET OF 4 
RINGS FROM THEIR SHEAR CENTER 4 

THE MOMENT OF INERTIA OF THE CROSS-SECTION OF THE 4 
K TH SET OF RINGS ABOUT THE Z-AXIS PASSING THROUGH 4 
THEIR CENTROIO 4 

THE MOMENT OF INERTIA OF THE CROSS-SECTION OF THE 4 
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♦ K TH SET OF RINGS ABOUT THE X-AXIS PASSING THROUGH ♦ 

♦ THEIR CENTROID ♦ 

♦ GJR ( K I THE TORSIONAL STIFFNESS OF THE K TH SET OF RINGS ♦ 

♦ ♦ 


♦ THE REMAINING PARAMETERS OF THIS PROGRAM ARE DESCRIBED IN THE ♦ 

♦ C IC T IONAR Y OF VARIABLES ♦ 

♦ ♦ 

* SUBROUTINES REQUIRED ♦ 


♦ 

INTGRL 

4 

♦ 

xx ! 

4 

♦ 

GAUSS 

4 

♦ 

SHELLl 

4 

4 

SHELL2 

4 

♦ 

R I NG 1 

4 

♦ 

PING2 

4 

♦ 

RING3 

4 

♦ 

RING* 

4 

4 

RINGS 

4 

♦ 

R I NG8 

4 

♦ 

EIGEN 

4 

♦ 

JACOBI 

4 

♦ 

MATMUL 

4 

4 


4 

4 

FUNCTION SUBROUTINES REQUIRED 

4 

4 

RSHL 

4 

4 

RRRT 

4 

4 

RSHLT 

4 

4 


4 


INT EGER NBC,BC(2,4), BUI 2). TITLEK 7),TI TLE2I 71 

DIMENSION T{ 1 ,161 ,PSI1» , ESI II ,ASU) »ZIS(1 ) .Z2SII I, VIS (ll,Y2S( 1 ), 

1Z IS <11, Y ISC I ),YZ IS(l),GJS(l),RX( 2,1 1» .PR I 2) ,ER( 2) ,AR( 21 ,E1R< 2) . 
2E2R(2I,ZIR(2I , X I R 1 2 » , GJ R 1 2 > ,X IS , 5 S ) , XXX < 2, 2, 55 ) , C I R ) . 

3 SUM! 181 ,NNLI1I,NNM2),ST<75),TSU.42I,SS(1,30> 

4,NNR(2).NR(2.2).CR(2,AO),RIC2, 5AI.RCGI2I 
DIMENSION AM 90,90) , AM (90 .90 I ,VECR(90,90 ),EVR(90 ),LCI90) 

DIMENSION XXXXl 8100) , Y( 81 00) , Z I 81 00) , MC I 90) , E VI I 90 ) ,1 NDICI90) 
COMMON DR (9 ) , RI9 ) , DRV I 5 ) , RV ( 5 ),R 1 C 8 ) , RR 1 1 8 ) , R 2( 1 0) ,RR 2( 1 0) ,R 3( 2 ) , 
1RR3I2) ,RA<5) ,RRA( 5 I ,R5( 181 ,RR5I18) , R6 1 1 1 ) ,RR6( 1 1 ) , P | , XK , A A, X l< 5 ) , 
2XRI 2 ), EIRK, E2RK,N,NB ,NBC * K » KB ,NSA 
CATA BC/IOHCLAMPED-FA.IOHEE , 10HF REE LT SUP , 1 OH PORTED ,10 

1HCLAMPED CL.IOHAMPEC , 10HFR EE-FRE E , 10h / 


EQUIVALENCE THE STIFFNESS AND EIGENVECTOR MATRICES 


EQU IVAL ENCEI AMI l, 1 1, XXXXl 1) ) 

E QUI VALENCE I AK 1 1,1 ),Y ID) 

EQUIVALENCE! I NOIC 1 1) , LC 1 1 ) I 

EQUIVALENCEIVECRI1, l)»Z! l),Xtl, t ) ) , I VECR I 1 , l 1 ) , X XXI 1 , 1 , 1 ) ) 

EXTERNAL SHELLl 

EXTERNAL SHELL 2 

EXTERNAL RING1 

EXTERNAL RING2 

EXTERNAL RING3 

EXTERNAL R INGA 

EXTERNAL RINGS 

EXTERNAL RING6 


KRRR SHOULD BE EQUAL TO THE FIRST DIMENSION OF AK, AM, VECR 
MATRIC1ES 


KRRR = SC 
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10000 2ERC=0.0 
NE X l T =0 

REAC THE NAME 0E THE B0UN0AR Y CONDITION 

REA0(5,8»BCR,NQUIT 
8 FORMAT! 2A10,59X,ll) 

f F (N0U1T .NE. 0) CALL EXIT 
2CCe WRI TEI6, 10011 

1001 FORMAT! 1H1.6X, 67! IH*»,//,6X, 65HFREE VIBRATIONAL ANALYSIS OF STIFFE 
INFO OR UNSTIFFENEO CIRCULAR OR, / , 1 3X , 51 HNONC IRCUL AR CYLINDERS WITH 

2 ARBITRARY END COND 1 T ION S , // » 6X , 67! IH* I , //// I 

IDENTIFICATION OF THE 80UN0ARY CONDITION AND ASSIGNING A CODE 
NUMBER NBC AS FOLLOWS 

NBC = l FCR CLAMPEC-FREE 

NBC = 2 FOR FREELY SUPPORTED 

NBC = 3 FOR CLAMPEC CL AMP EO 

NBC « 4 FCR FREE-FREE 

CO 2 J=l,4 
DO 3 1=1,2 

IFIBCRIII .EQ. BCII.jn GO TO 3 
CO TC 2 

3 CONTINUE 
NBC= J 

GC TC 4 
2 CONTINUE 

WR I T E ( 6, 2003 I 8CR 

2 0C3 FORMAT ( / / , 1 X ,19 H***** ERROR ***** , 2 7HB0UNDAR Y CONDITION REAO IS 

1 ,2A10,/,20X ,46HTHt BOUNDARY CONDITION MAY NOT bt wORDEi) RIGHT, /20 
2.37HOR THIS BOUNOARY CONDITION MAY NOT B E , / , 20X , 2 5HA VA IL ABLE IN TH 

3 IS PROGRAM ) 

NE X I T = 1 

READ ANO WRITE THE GENERAL INFORMATION 

4 PEACI5.60 )NG,KG,LL,NL,KK,NK,MMIN,MMAX,MSA,NM!N,NMAX,NSA,NEO,IR, 

1 NWK , NWM ,NWEV 

fcC FORMAT! 2014) 

WR IT E (6, 10015 ) NG,KG,LL ,NL,KK,NK,MMIN,MMAX,MSA,NMIN,NMAX ,NSA ,NEO, 

1 I R , NW K , NW M , NW EV 

10015 FORMAT! 26X, 25HGENERAL INPUT l NFORMA T l ON , / , 26X ,25 ! 1H- ) , / , 8 X , 

16HNG =,I4,2X,6HKG = ,I4,2X,6HLL =,I4,2X,5HNL =,I4,3X, 

26HKK =,I4,/,8X,6HNK = , 1 4 , 2X ,6 HMM l N = , 14 , 2X , 6MMM AX =,I4,2X, 

35HMSA =, 14, 3X,6HNMIN = , I 4 ,/ , 8 X, 6HNMA X = , I 4 ,2 X *6HNS A =,I4,2X, 
46HNE0 =,I4,2X,5HIR = , 14, 3X , 6HNWK = , I 4 , / , 8 X, 6MNWM =,!4,2X, 

56HNWEV =,I4, ////) 

PI=3. 141502653585703 
F.I2=PI*PI 

174 IF! NL .GT. LL) GO TO 176 
GO TO 177 

176 WRI TE !6 »l 781 NL.LL 

178 FORMAT! //, IX, 10H***** ERROR ***** , 5HNL = ,I4,5X,5HLL = ,14,/, 

1 20X,28HNL CANNOT BE GREATER THAN LL » 

NEXI T= 1 

177 IF ( NK ,GT. KK I GO TO 170 
GC TC 180 

175 WRI TEI 6*181)NK,KK 

181 FORMAT!//, IX, 10H***** ERROR ***** , 5HNK = ,I4,5X,5HKK = ,14,/, 

l 20X ,2 8HNK CANNOT eE GREATER THAN KK I 


14 



non <~>n o n o n o ooo 


NEXIT- l 

COMPUTE THE CRDER OP THE MASS ANO ST 1PPMESS MATRICES 

180 IF1NEXIT .GT. 01 GO TO 18000 
MD=1 

1FIMSA .LT. 21 MO-2 
MS= ( MMAX-MM1N )/ M0*1 
ND=1 

I F ! N EO .LT. 21 NO-2 
NS=!NMAX-NMIM/NO*l 
MN=M S*N S 
10=0 

IF! NBC .EQ. A .AND. MSA .NE. 1) IO-NS 
IFINMIN .GT. 01 GO TO 20A5 
IFINSA) 2046.20A6.20A7 

2046 MN3 =3*MN-I O-MS 
GO TO 2048 

2047 MN3=3*MN-I0-2*MS 
GO TO 2048 

2045 MN3= 3*MN- 10 

ZERO OUT THE UPPER TRt ANGULAR MATRIX OF MASS AND STIFFNESS 
MATRICIES 

2048 DO 2C04 I=l,MN3 
CO 2004 J= I »MN3 
AKI I ,J1=0.0 

20C4 AMI I,J)=0.0 

READ ANC WRITE THE SHELL DATA 

P EA C ! 5. 1000 1 T ITLE 1* T I TLE 2 
1009 FORMAT !7A10«/ • 7A 10 ) 

WRITE(6,1003ITITLE1,TI TLE2 

1003 FORMAT ( 29X , 19HS HELL OAT A , / , 29X, 191 IH- ) . // // . 5X ,7 A10 ,// .5X 
1 .7 A10 »////) 

READ!5,65IPC,EC,XNU,H,AA 
65 FORMAT (5E15 .81 

WRITE! 6, 1002 1 PC ,EC , XNU,H ,AA ,BCR 

1002 FORMAT! 10X, 12HMASS DENS IT V, 10X, 2H= .E15.8.18H LB SEC . **2/I N. ** A / / , 
110X.2AHM00ULUS OF ELASTICITY = , E15 .8 , 10 H L B/ IN . **2// , 10X, 1 5HP0 I SS 
ZON'S RATIO, 7X.2H- , E 1 5. 8 , // , l CX.9HTH I CKNE SS , 1 3X ,2H = .E15.8.7H INCH 
3ES,//, 10X.6HLENGTH, 16X, 2H= .E15.8.7H INCHE S, //, 10X, 14HEND CONDITIO 
A NS » 8X . 3H= . 2 A10 ) 

PC=PC*H*2.0 
I P ( LL .EO. 0) GO TO 85 

PEAO ANO WRITE THE STRINGER DATA 

WRI TE < 6 »l OOA l LL.NL 

100A FORMAT! 1H1, 26X, 25HS TRINGER OAT A , / .27X.25 1 1H-I ,// ,17X , 
l A3H ( THE UNITS ARE SAME AS THOSE OF SHELL DATA I.//.23X .2 8H TOTAL NUM 
2BER OF STRINGERS = .IA./.15X.A1 HNUHBE R OF DIFFERENT KINDS OF STR IN 
3CERS = , IA./.5X.67I 1H-I1 
IZ1 =0 
IZ 2=C 
IYI=0 
I Y2 =0 

DO 66 L= 1 »NL 
R EAC ( 5. 60 1NNL (L I 
NNM=NNL!L) 
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RE AO (5, 65HT(L,I ), I* l, NNNl I 

PEAO (5,65 IPS (LI. ES(L I , AS II I . Z IS (l I, Z 2SIL I , VIS! L I , Y2S( L I , Z I S( L > , 
1YIS(L) ,YZIS(L) ,GJS(ll 


IFI21SIL I 

.NE. 0.0 

1 

121-1 

IFIZ2SIL) 

.NE. 0.0 

) 

I Z2-1 

IF(YISIL) 

.NE. 0.0 

1 

I Yl«l 

IF !Y2S (L I 

.NE. 0.0 

) 

IY 2* 1 


WRITE 16 ,1005) NNKLI ,PS!L1 ,ES III ,AS(L),Z1S(L I, Y1S (L > , Z2S (l I, Y2S( L I , 
1ZIS(L>,YIS(L),VZIS(1>,GJS(L) 

1005 FORMAT!//, 15X, 14, 41H STRINGERS WITH THE FOLLOWING PROPERTIES ,//, 
15X.18HMASS OENSITY ■ ,E15.8,2X,17HW0D. OF El AS . = ,E15.8,/,5X, 

24HAREA, 12X, 2H= , E 15 . 8 , 2 X , 17HSHE AR CTR. 1211* f El 5. 8, / ,5 X , 1 8HSHE AR 
2CTR. (Yll = ,E15 .8,2X,17HCFNTP0I0 IZ2I « , E15.8, /, 5X, 18HCENTR01D 
4 I Y2 1 = ,E15.8,2X,17HIN£RTI A UZZI * , El 5 . 8 ,/ ,5X , 1 8 HINERT I A ( IY 

5Y ) = ,E15.8,2X, 17HPR00. INER .1 IYZ »= ,E 15. 8, / ,20X, 22HTORSI ONAl ST I F 

6FNESS = ,E15.8,//,5X,43HL0CATE0 AT FOLLOWING THETA VALUES (DEGREES 
71,/ ) 

WRITE(6,1006)(T(L, II, I* 1 , NNNl I 
lOCfc FORMAT! 4X, El 5. 8, IX, El 5. 8 ,IX ,E15 . 8 ,1 X , E15 .8 1 
CO 2000 1=1, NNNL 
T(l,l)=T!L»l)*PI/0.18 E*0 3 
20CC CONTINUE 

WRIT E (6 , 1010 ) 

1010 FORMAT !/,5X«67(lH=ll 

COMPUTE THE MOMENT OF INERTIAS WITH RESPECT TO AXES PASSING 
THROUGH THE SHEAR CENTER OF STRINGERS 

IFIZ2SILI .EQ. 0.0 I GO TO 182 
ZlS(LI=ZIS(l)>AS(l)*Y2S !LI*Y2SIL) 

182 IFIY2SIL) .EQ. 0.0 I GO TO 66 
YlS(l>=YIS(L)*AS(D*Z2S(l l*Z2Sll I 
YZISILl=YZIS(L»»ASm*Y2Sm*Z2SIU 
66 CONTINUE 

REAO AND WRITE THE RING CAT A 

85 IFIKK .EQ. 0 » GO TO 86 
WRI TE(6 ,10071 KK.NK 
1007 FORMAT! 1HI, 30X.17HR I N G OAT A , / ,31 X ,1 7 ! 1H-I ,// , 17X ,43H (THE U 
1 NITS ARE SAME AS THOSE OF SHELL OAT A),//,24X, 24H TOTAL NUMBER OF RI 
2NGS = ,I4,/,17X,37HNUMBER OF DIFFERENT KINDS OF RINGS = , I4,/,5X, 
367! 1 ) I 
I El -0 
IE 2*0 
NKR= 1 

00 75 K=1 , NK 
READ! 5, 60INNKIK) 

NNNK=NNK (K I 

PEAO (5,651 (RXIK.I I , 1=1 , NNNK I 

READ! 5, 65IPR(K),ER(KI,AR(K> ,E 1R ( K » , E2R ( K t ,ZI R(KI , XI R(K) , GJR(K) 
IFIEIR(K) .NE. 0.0 I I E 1* 1 
IF(E2R(K» .NE. 0.0 I IE2-1 
CG=EIR(K)*E2R(K) 

IF ! K .EC. II RCG!L)*E1R(1)+E2R(1I 
00 10008 1=1, NKR 
IF(RCGII) .EO. CGI GO TO 10009 
10008 CONTINUE 
NKR=NKR*l 
RCGINKR l = CG 
100C9 CONTI NUF 

WRt TE( 6 , 1008) NNX! Kl ,PR(KI ,ER(K) , AR! K) ,EIR(K),E2R(KI,ZIR(K),XIR(K), 
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1GJRIKI 

10C8 FORMAT (//.17X.I4.36H RINGS WITH THE FOLLCWtNG RROPERT I ES, //« 5X« 18H 
1MASS DENSITY • ,E15.8,2X,17HM0O. OF EUST. * .615. B/.5X.4HAREA, 
2 12X.2H- .E15.8f2X.l7HS HEAR CTR. (Ell» , E 15. 8/, 5X, 18HCENTROIO 1621 

3 - , E15. 8. 2X.17H INERTIA (1221 » , E15.8/ ,5X .18HINERT IA ( IXX ) 

3.E15.8.2X. 17HTORS. STIF.(GJ)- .E15.8.//.5X, 38HLOCATEO AT FOLLOWING 

4 X VALUES (INCHES I ,/l 

WRITE(6,1006MRX(K,I1 ,1-l.NNNKI 
WRIT E (6, 1 01 0 ) . 

IF(E2R(K) .60. 0.0 I GO TO 75 
X(R(K I* X IR (K I *AR (K l*E2R(Kl*E2R(K) 

75 CONTINUE 

CENTROI OAL INFORMATION OF RINGS 

00 1C010 I *1 »NKR 
NNR ( I 1=0 
CC 10011 K= 1 » NK 
CG=6 IR(KI ♦E2R( K) 

IF( CG .NE. RCG(II) GO TO 10011 
NNR ( I ) =NNR( I ) *1 
NR ( I ,NNR( IM»K 
0011 CONTINUE 
0010 CONTINUE 

86 IFIMMIN .GT. 01 GO TO 67 
I F ( NBC .LT. 4) GO TO 99 

IF MMIN = 0 INCREASE MMIN AND MMAX BY l 

MPIN=MMIK*1 

mmax=mmaxm 

67 I MO 
NCHNG-0 

IFINMIN .GT. 01 GO TO 2040 
NMIN=1 
NMA X=NMA X*1 
NCHNG= 1 
2040 CONTINUE 

EVALUATE THE LONGITUDINAL INTEGRALS AND THE X( OUTPUT OF 
SUBROUTINE XX » VALUES ANO STORE THEM 

CO 70 M=MMIN,MMAX,MD 
K = P 

IF ( M 5A .NE. 1 .AND. NBC .60. 4) K=M-1 

CO 70 MB-M, MMAX.MD 

KB«PB 

IFIMSA .NE. 1 .ANO. NBC .EO. 41 KB-MB-1 
IM— IMM 
CALL INTGRL 
CO 80 1=1,5 

x(i.iMi=xim 
8C CONTINUE 

IFIKK .EO. 0) GO to 70 

00 71 I =1 .NX 

XXX( l.I ,IM»=0.0 

XXX ( 2, I, IM 1=0 .0 

NNNK =NNK 1 I ) 

00 71 KKKK* l .NNNK 
XMRXI I.KKKK > 

CALL XX 

XXX( 1 . 1 , I Ml =XXX( l.I ,1 Ml ♦ XR( 1 1 
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XXX (2. I , 1 N )*XXX 1 2 « I . INI tXRI 2 ) 

71 CONTINUE 
70 CONTINUE 

evaluate the constants of the shell equations 

C*EC*H*H*H/(12.0 *(1.0 -XNU*XNU)» 

S5=2.0*D 

S7*S5*XNU 

S3 *0*1 1.0 -XNUI 

S6=3.0E0*S3 

S8=4 .0E0*S3 

S1=12.0EC*S5/IH*H) 

S4=S 1*XNU 

S2=S8*3.0E0/ IH*HI 

IF ( LL .EQ. 01 GO TO 167 

EVALUATE THE CONSTANTS OF THE STRINGER EQUATIONS 
00 35 1*1, NL 

tsii,ii*psii>*asiii 
T s<i,2»*psm*zism 
TS( 1 , 3 l*PS( I ) *Y I SI I » 

TSI I,4)*TSI 1,2 >»TSI 1,31 
TSI I , 5) = TSI I , A) ♦ TSI I ,41 
SSI I,1) = ES( 1 1 *A S < 1> 

SSU .21 =ESI I )*Z IS 1 1 ) 

ssi i , 3) =e si n*Yi sin 

IF (USUI .EQ. 0.0 I GO TO 15 
TSI I ,61 «2.0E0*TSII,1)*Z1S(I> 

TSI I«10)*TSII,2)*Z1SII) 

TSII,7I*TSI I,10)*Z1SI I) 

TSI I ,BI *TS( I ,10) ♦TSI 1,10) 

TSI I,9)*TSI I,1)*Z1SI I )*Z IS! I ) 

TSI I ,11 ) =TS I 1,7) *TS 11,7 ) 

TSI I,12) = TSI I,9)*TSII ,9) 

SSI I « 4 ) = SS I I, 1 ) *Z 1SI I ) 

SSI I ,7 ) =SS 1 1 ,2 ) *Z1S I I ) 

SSI I , 6 ) = S SI I , 7) ♦ SSI I ,7) 

SSII»5)=SSlt»7)*ZlSI I) 

SSI I ,8) *SSl I ,4 ) *Z1 SI I) 

15 IFIZ2SII) .EQ. 0.0 ) GO TO 20 

TSI I »13)=2.0E0*T$I I, l )*Z2S 1 1 ) 

TSI I ,14) = TSI I ,13)*ZISII ) 

TSI 1,15 ) = TS I I, 14 1+TSI I, 14) 

SSI I ,9) =SSI! ,1)*22SII) 

SSI I »10)=SSI I »9)*Z1SII)*2.0 
20 IF IY IS I I ) .EQ. 0.0 ) GO TO 25 

TSI l ,281 =TS 1 1 ,1)*Y1SII) 

TSI I, 16 >*TSI I, 2 8>* TSI 1,28) 

TSI 1,17 !*TS I 1,16 >*Z2 SI I > 

TSI I ,181 = TSI 1 *28 1* Y1 SI I ) 

TSI I, 251* TSI 1,18 ) ♦ TSI 1,18) 

TSII,26)=TSII,3)*YlSm 
TSI I ,19I*TSI I ,26)*V1SII I 
TSI 1,24 )*TSI 1,19) ♦TSI I, 19) 

TSI I ,20) =TS 1 1 ,25 ) *Z2 S 1 1 ) 

TSI I,2i)=TSI l, 16>*Z 1SU) 

TSI I,22)*TSII,26 I ♦T S I 1,26) 

TSI ! ,23)*TS(I ,21 )*Z2 SI I ) 

TSI 1,27 )*TSI I,23)/2.0 
SSI I,ll ) *SS 1 1 , 1 ) *Y1 S I I) 
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SSII,12)-SS< I, ll)*Z2Sll 1 
SSI 1,13 l-SS 11,11 I*Y1SIII 
SSI l ,17)-SS< I ,3I*YISI I) 
ssi i,i4i*sst t.mtvisu ) 
ssi i,i4)*ssti,ii )*zisii) 

SSI 1,18) -SSI I , 121*21 SI I ) 

SSI 1,19 )-SSt I,12)*Y1SII ) 

SSI I ,15) >SSH,l9)*2.0 
25 !F| Y2SI f I .ED. 0.0 ) GO TO 30 

TSU,42)-TSIl,l)*Y2St I) 

TSI I ,29)-TStI,42)*TSII,42) 

TSI 1,301-TSI t,29l*ZlS<tl 
TSII,31)-TSI!,2«»>*YISI!) 

TSI I *321 •TSI 1 ,311*21 St I ) 

TSI 1,35 I- TSI I,30)*ZISI! ) 

TSII ,39)«TSII,31)*TSII,3l) 

TSI 1,40)-TSI 1,351/2.0 
TS{!,36)*2.0*PStt)*YZIStI) 

TSI1 ,33) -TSt! ,36)*Y1SII> 

TSI l,34l«TSt 1 , 331*2 1 SI I > 

TSI I #37 )*TS 11,36 )*ZISI 1 1 
TSI I , 381-2. 0*TSt I ,34 ) 

TSI I»41)*PSI I »*Y2 1 SI 1 1*2 1 SI 1 1 
SSI (,20)»SS(!,l>*Y2Sin 
SSI I ,21 1 -SSI 1 ,20)*Z1SII I 
SSI I, 221* SSI I,20)*2.0E0*YlSf I) 

SSI 1, 23) -SSI !,22>*Z1SIII 
SSI 1,241-SSIf ,23)/2. 0 
SSI 1,25 )*SS II,21)*ZlSin 
30 IFIVZISIl) • EQ. 0.0 ) GO TO 35 

SSI 1,281-ESI I )*Y2ISI I) 

SSI 1,29 l-SSI 1,28 I*21SI II 
SSI I ,30)* SSI T,28)*Y1SII I 
SSI 1,26 l*SSt I • 30 )+SS I 1,30) 

SSI I ,27 )*SSI I ,26 )*21S (I ) 

35 CONTINUE 

167 IFIKK .60 . 0) GO TO 168 

EVALUATE THE CONSTANTS OF THE RING EQUATIONS 

00 125 K-l,NK 
CR(K,1)*2.0E 0*6 R I K) *ZIR I Kl 
CRIK,Z)-2 .0E0*ERIK ) *AR I K ) 

CRIK, 3) *2.0E0*ER IK )*X IRIK) 

CRI K , 21 )*2. 06 0*G JR I K) 

CR IK ,22 )*2.0E0*PRIK l*AR IK 1 
CRIK, 23) -2.0EO*PR(K)*ZIRIK) 

CRIK,24)*2. 06 0*PRIK)*XIRIX) 

IFIE1RIK) .EQ. 0.0 ) GO TO 126 

CRIK ,4) «CRIK,1)*EIRIK) 

CRI K » 9 )*CRI K * 21*6 IRIK) 

CRtK,5l*CRIK,9)*ElR{K) 

CRIK,6)-CRIK,3I*E1RIK)*E1RIK) 

CRIK.7 l-CRI K, 9 ) ♦CRI K ,9) 
CR|K,8)*2.0E0*E1RIK)*CRIK,3) 

CRIK, 10)*CR|K,4I*E1RIK) 

CR|K,25)*CRIK,21 )*E1RIK) 

CRIK, 26) *CR| K, 25) *6 IRIK) 
CRIK,ZTI«CR|K,25)4CRIK,25) 

CRIK, 31 )*CR|K,22 )*61 RIK ) 
CRIK,30)*CR|K,31)*E1RIK) 

CRIK, 28 l*CRIK,31 1+CRtK, 31) 
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CRI K ,291 =CRIK,23 »*2.0E0*E1RIK> 

CR CK, 32 »«CRIK,24I*E1RIK >*E1RIKI 
CRIK,33I=CRIK,24 )*2 .0E0*EIRIK» 
CR(K,34»=CR!K,23»*E1RTK»*EIRIK» 

126 IF ( E2R (K I .EQ. 0.0 ) GO TO 127 

CRIK,14I*CRIK,2>*E2RIK) 

CRIK.ll »=CRIK,14I*CR(K,14I 
CRIK,16)=CRIK,14)*E1RIK I 
CROC. 17I=CR1K,16UCRIK, 161 
CRIK, 12>=CRIK, 17I*E1RIK» 

CRIK,13) = CRIK,17 l*CRIK, 17 I 

CR1 K ,151 =CRI K,17|»CRIK,16I 
CR IK, 36 I=CRIK,22»*E2RIK> 

CR(K«35)=CR(K.36)+CRIK,36) 

CR(K,40>=COIK,35)*E1ROO 
CRIK.38»“CR(K,40»*CRIK, 401 
C R ( K ,391 =CR I K ,40 1 *CR I K ,38 I 
CR I K , 37 1 =CR I K , 40 1 *E 1R I K ) 

127 IF ( I P .EQ. 01 00 TO 125 
CRIK,19»=CRIK,3)*EIRIK» 

CP(K,18»=CRIK,6)*CRIK,6» 

CRIK ,20 ) =CR IK,16 >*E1RIK » 

125 CONTINUE 
168 N0C x 0 

THF 00 LOOP OF N 

DC 90 NASA=NMIN,NMAX ,N0 
N=NA SA 

IF ( NCFNG .N6. 0 IN=NASA- 1 

A0C*N0CM 

AN* FLOATINI 

AA2* AN* AN 

AFC = C 

THF CC L COP OF N 8 

DO 91 N A SB = NX I N,NXAX,ND 
AB= A AS 8 

1FINCHNG .NE. 0>NB=NASB-1 

ON* FLOAT! NB I 

AEC=NECa1 

BN2=8N*8N 

ABN= AN* BN 

ABN2 =ABA* ABN 

ARNB=ABN*BN 

ABNA=ABN*AN 

IFILL .EC. 0) GO TO 169 

EVALUATE THE CIRCUMFERENTIAL OUANTITIES CF THE STRINGER 
75 IS THE TOTAL NUMBER OF TERMS IN THE STRINGER ENERGIES 

DO 94 IL* 1* 75 
ST ( I L 1*0 .0 
S4 CONTINUE 

CO 95 L* l.NL 
AKMNNUL) 

00 l Cl Kl=l,8 
1C1 C. IK l 1=0 .0 

CO 96 L 1=1, NAN 
TN=AN*T! L ,LI I 
TNB=BN*TIL,L IV 



I F ( N S A .EO. 1 » GC TC 97 
CN= COSCTN) 

CNB = COS CTNBt 
SN = SIMTNt 
SN8= SINCTNB) 

CC TC 98 
S 7 CN= St N( TNI 
CNR= SINCTNRt 
SN - COS(TN) 

SNB= COS( TNB) 

98 CC=CN*CNB 
CS=CN*SNB 
SSS=SN* $N0 
SC= SN *CNR 
SR=PSHL <T CL ,1 Ctl 
SR2=SR* SR 
CCI)=CC11»CC 
C <3 1 =C C 3 I *SSS 

IFC21SCL1 .FO. 0.0 » GO TO 96 

C (2 1 = CC 2 UCS 
CC4I *CUI»SC 
C C 5 1 =C C St ♦ SSS/SR 
C C6 J =C C 6 1 -»SS S / S R 2 
CC 7t =CC 7UC S/SR 
C C 8 1=CC 8 ) ♦SC /SR 
96 CONTINUE 

ST( 1 1 =ST< l)*SSIL,n»CIII 
STC2l=STC2)*SSCL„2!*CC3t 
STC3t=ST(3t*SSCL«3i*CCl ) 

ST( 7 3 1 = S T ( 73UGJSa.) *ABN*CC6> 

ST {74 l=ST < 74 t ♦GJSit 1 *CC 6) 

ST I 75l*STC75)*GJS(Lt*BN*CC6t 
srt 2<5 li=sr I 29?) <-TSC L 0 1 » *C I 1 ) 
ST(30)=ST{30}vTS8L,2)*C!3> 

STI 31t = S70U*T$«L,lt*C I 3t *TS C L ,4 1* C C 6 > 

ST< 32 }= ST (32 J<vTS<L,5)*BN*CC 6) 
ST(33t=STC?3)4TSCL*3l*CCl ) 

ST< 34 1= ST J 34t»ABN*TSCL,4t*C(6l 
ST(35t = ST(38)*TSU, l )*C< 1 t 
IF ( 21 SC It .EC. 0.0 > GC TO 102 

STC AUSTCAt-SSCL, AM>CC 1 ) 

ST (5 t = STC5 t ♦SS(L P 5>«CC6 l»SSIL,6 t*C( St 
ST( 61 =STC6t*BiM*< S S 1 1 ,51 *C (6 > +SS C l ,7 ) *C C 5 I » 
ST(7J=ST(7t*SS(Lf8t*CC l) 

ST(8 t =ST(8 ! »SS( l ,5 »*C(6 t*ABN 
S T ( 3 6 I = S T { 36)-TSCl,6l*CClt 
ST< ?7t = $T{ 37 toTSIL ,7t*C< 6 H-TSCL , 8»*CI 5t 
ST ( 38 t =STC38l»TSCL,9t*C C6UTS CL ,6 )*CC5 t 
ST(19)=STj 39 1 »C TS(L»lll*C(6l*TSCL»81*C<5tl*8N 
ST<40l = ST!4Pl-MTS<L,l2)*Cl6im<L,6>*CC5t ) *BN 
ST( All =ST5 41 f » TSC L ,9»*C Cl t 
ST <42 )=ST? 4 2 H4BN*TSCL, 7|*C< 61 
STC43 l=ST(435«49N»TSCL,91*CC6l 
1 C2 IF C 2 2SC L I .EO. 0.0 I GO TO 103 
STC9l = ST<9t-SS<U9!*C(l) 
STtl0t=STC10t*SSCS.,10t*CCl! 

STC 44) = STC 441-TSCL, 13)*CC 11 
ST C45t = ST C451ATSCL, 131*CC51«-TSCL, 141 *CC 6 1 
ST<46l = STC46l* C 7SCL,13)*CC5 MTSC L,15)«Ct6)UBN 
STC47>=STC47)*A3NfTSCL, 14UCC6I 
STt481=ST(4Bt*TS«C,14l*CCll 
1C3 I F C Vi SC L 5 .EQ. 0.0 ) GO TO 104 
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STIin=STIll»-SS(L, 11»*C(2»*SS(L,12I*C« r» 

ST( 12 )= ST| 12 ) ♦BN* SSI L ,1 2 ) *C ( 7) 

ST(13)=ST( 1 3 * *C I 3)*SSCL, 14>*C( 6)^SS(L, 151*01 51 

ST{ 14 ) =ST I14|* (SS (L, 12) »SS (L ,16 I )*CI4 l-( SSIL, 17 |>SS(L #181 )*C( 81 

ST( 1 5 » = ST* 15)*BN*(SS<L,14)*C(6)-SSCL,19»*CC5)) 

ST(16)®ST(16) 4ABN*SS IL » 14 )*C( 61 

STI 17)=ST(l7)-< SSI L» 17) +SS (L , 1 3 » ) *1 BN*Cf 7 ) »AN*C 18) > 
ST(49)=ST(49)-TS(L, 1 6I*C ( 2 ) ♦ TSI l, 17 ) *C ( T ) 

ST ( 50 ) = ST (50 )*BN*TS(L«17 ) *C ( 7 ) 

ST(51 J=STC5l)*TSCL,18>*C<3UTS(L,l9)*CI6l-TSIl,20)*C(5) 

ST ( 52)=ST(52 I ♦TS( L • 1 8) *C I 6) 

Sr(53)=ST<53)4<TS<L,2l)*TS(L,l7l>*C<4)-(TS(L,22t*TS(l,2m*C(8) 

ST* 54)*ST(54»*( TS(lt24)*C(6l-TS(l,20)*C(5l)*BN 
ST(55)=ST(55I*TS(L,25)*BN*C(6) 

ST156I *ST(56I-TS(L,16)*C(8) 

ST( 57)=ST( 57)*ABN*TSIL, 19)*C(6> 

ST (50 I = ST (5 8 )-(TS(L, 26 )*TS( L.27H *(BN*C( 7) »AN*C( 81 ) 

ST ( 59l=ST(59)*ABN*TS(L,18)*C(fcl 

ST (60)= ST (60 )-TS(L ,28)*<BN*C( 7)*AN*CI 8) ) 

1C4 IF (Y2 S ( L ) .EC. 0.0 ) GO TO 105 

ST( 18)=ST( 18)-SS(L,20)*C<?)-SS(L,21)*C(7) 

ST( 10 1 = ST( 19 )-8N*SS(L, 2 1)*C( 7) 

ST(20)=ST(20)*SS(L»22)*C<3I»SS<L,23 )*C(S ) 

ST( ?1)=ST( 21)*SS(L,24)*BN*C( 5) 

ST(22)=ST(?2)*SS(L,2l )*C ( 4 ) ♦SSC l , 25 1 *C( 8) 
ST(23)=ST(731*SS(l,25)*IBN*C(T)*AN*CI811 
ST(6l )=ST( 61 l-TSIL ,29)*C( 2) - TS( l , 30 1 *C ( 71 
ST (62 1=ST(62 >-BN*T S ( L , 30 ) *C ( 7 ) 

ST( 63) = ST( 63) *TS( L« 31 ) *C ( 3) ♦ ( TS ( L , 32 ) - TS ( L , 3 3 ) ) *C (5 I -TS f 1, 34 |*C('61 
ST (64 ) = ST (64)tTS(L. 31 ) *C ( 6 ) 

ST( 65) =ST (65) «• < TS(L,30I ♦TSIL »36 ) )*C(4 >♦( TS <L, 35 )*TS (L,37) )*CI8) 
ST(66)=ST(66)-BN*( ( - T S( L , 32 ) ♦ TS( L ,33 ) ) *C ( 5) * TS ( l ,30 ) *C (6 ) I 
ST { 67 ) = ST (67 )-TS(L«29)*C(8 ) 

ST( 68)=ST(60)*BN*TS( L,39I*C(6) 

ST ( 69 )= ST( 69 ) ♦( TS(L*40)*TSI L»4l) ) * ( BN*C ( 7) «-AN*C(fll ) 
ST(70)=ST(70)-ABN*TS (1,34 )*C(6) 

ST( 7 l ) * ST( T1)*ABN*TS(L,31I*C(6) 

ST (72 )=ST<72)-TS(L,42)*(BN*C(7) ♦AN*Ct 8) ) 

105 IF(YZIS(L) .EQ. 0.0 ) GO TO 95 

ST( 24 ) = ST( 241-SSIL ,26)*C( 5) - S St l , 27 » *C 1 6 ) 
ST(25)=ST(25)4SS(L,28)*C(4)4SS(L,29)*C( 8) 

ST(26) = ST(26)-BN*( SS(L,30)*C(5)*SS( L,2T)*C(6) ) 

ST( 2 7)=ST(27 (♦SSIl. ,29)*(BN*C( 7)*AN*C( 8) ) 

ST(28)=ST(28)-ABN*SS (L,27 )*C(fcl 
55 CONTINUE 

INTEGRALS OF SHELL 

169 CALL GAUSS ( NG ,KC. ,7ER0 , P I , 9 ,PHI ,0R , SUN ,R , SHE L LI ) 

00 131 10=1,9 

I E ( ABS(RIIO) ) .LF. 0.1F-08) R(IC>=0.0 

131 CONTINUE 

IF{ IP ,?C. 0) GO TO 130 

call GAUSS INuiMi, L tKo , P l, 5, PHI ,Ui<V, SOM,n V, iHtLLC) 

OU 132 lw=l ,5 

I F < ABS(RVIIQ)) .LE. 0.1E-08) RV(IQ)=0.0 

132 CONTINUE 

INTEGRALS CF RING 

13C IF(KK .EO. 0) GO TO 133 
OO 146 LA =1 ,NK 
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CO 146 LB 3 It 54 
14o hi I L A , L 6 )=0 «0 

00 1 3 o KKI 3 1 ,N KR 
NRKK 1 1 3 NK IKK 1 1 1) 

El R K = E 1 K I NRKK II ) 

E2KK 3 E2K( NKKkI U 
Kd 3 NNKlKKll 

CALL GAUiSING.KG.ZEKG . P 1 , 6 . PHI , K 1 , SUM , a* l ,R 1NG 1) 
00 135 I0«1 tb 

IF I AaSIKR II 10 ) I .td. O.ie-Oo) RH1!IU) = 0.0 
CO l35 KR»1 , Kd 
NKKIKH=NM KKI ,KR) 

RIINKK 1KR, lOI^RKlI Id) 

135 CONTINUE 

IFIE1KK .Ed. 0.0 t GC TC 137 

CALL GAUSS I NG,aG« 2 ERO , P I , 10, PH l , * 2 , SUM ,RK 2 ,R INGZ I 
00 lib Id=l tlO 

IF C Ab SI HR 21 I U) 1 .LE. O.lt-Od! KR2IIG)=0.0 
CC 136 K R 3 1 * K d 
NKKi KR 3 NR(KKl , K k ) 

RIINKK IKK t Id *8 ) = Kk2 1 Id) 

136 CUNT INUE 

137 IFIE2KK .Ed. 0.0 I GC TO 134 

CALL GAUSS(No,KG,2cR0 , P 1 . 2t PHI , k i, S um ,kk i ,R I NO J I 
I F ( AbS ( kRj ( 111 ,Lt. 0 . IE-06 I RK3U) = C.O 
IF ( Ao SI HR i( 2) I .LE. 0.1E-06) KR3<2)=C.O 
CO 145 KR=l,Kd 
NKKI KK=NR(KKI ,Krl) 

K 1 ( NKK I KR till =K K 3 1 1 ) 

R I ( NKK IKK, 20 )=RR 3 1 2 ) 

145 CONTINUE 

134 IF ( IR .Ed. 01 GO TO lid 

CALL GAUSS (NGtK Gt ZERu ,P l,5tPHl ,R4, SUM,Rh<.,aING4» 
00 140 I 0=1 .5 

IF ( ABS(KK4(ld>> .LE. O.IE-Obi kK 41 ldl 3 C.O 
CU 140 KK 3 1 t KC 
NKKI KR 3 NK 1 KKI ,KKj 
R 1 ( NKK IK R i Id *20 )=Ka 4l Id) 

140 CCNTINLE 

IFIEIRK .Ed. 0.0 ) uL TC 141 

CALL GAUSS I NGtK Gt ZERO , P I , 18, Prt 1 1 R 5 t SUM ,R k 5 .R I NG5 ) 
00 142 I w 3 l ,16 

IF I A6 SI RR5I Id) ) .LE. 0.1E-06) HR5lld) 3 0.0 
OC l«t2 KR= 1 » k d 
NKKI KR =NR I KK I , KKI 
RI(NKKlKR,Id+25 ) =Rk5I 1m) 

142 CUNT INUE 

141 1FIE2RK .Ed. 0.0 ) GC TO lib 

CALL GAUSS! NG.KGtZERU , PI , 1 1 ,Prt I ,Ro , SUM , RRo , R1 NG6) 
CC 143 IU 3 ltll 

IF I ABSIRRbtldl) .LE. O.lE-Ob) RKb!Id)=0.0 

00 143 KR= 1 1 NO 
NKK IKH=NR (KKI , KR ) 

Kl I NKK IKK 1 1 0+43 ) 3 Rao I 1 w) 

143 CONTINUE 
139 CONTINUE 
136 CONTINUE 
133 18C=0 

C 

C Tt-E OU LOOP OF M 

C 

.00 S2 M = MMI N f UMAX , Mu 
IF I NChNG .Ed. 0) Gu Tu lOiO 
IF I NS A) 2031 ,2031 ,20.12 
2031 l 3 NDC«-lbC*NS-10 
. IEC=16C*1 

1 N-l+PN-I bC 
InN= 1 ♦MN*MN— H S 
GC TC 2033 
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203? I TEMP*N0C*I 8C*NS-I0 
I6C= I BC ♦ 1 
I = ITE*<P-I8C 
IN= I TEMP+MN-MS 
!NN= IN*MN-IBC 
GO TC 2033 

2C3C I =NDC* l BC*NS- l 0 
IBC= IBC *1 
IN® I +MN 
INN® I N*MN 

2033 JBC=0 

THF 00 LOOP OF MB 

00 <33 MB=MMIN,MMAX,MC 
IF ( NCHNG .EQ. 03 GO TO 2034 
IFINSAI2035, 2035, 2036 
2C25 J=NFC*JBC*NS-IO 
JRC=JBC*1 
JN=J*MN-JBC 
JNN=J*MN*MN-MS 
GO TO 2037 

2036 JTFMP=NEC*J8C*NS-IC 
JBC = JBC ♦ 1 
J = JT EMp- JBC 
JN=JTFMP*MN-MS 
JNN=JN*MN- JBC 
GO TO 2037 

2034 J=NEC* JBC*NS-IO 
JRC® JBC* 1 
JN=J*mn 

JNN= JN*MN 

2C37 IM=( IBC- l)*MS-IBC*I IBC-1I /2*JBC 
IFIjeC .GE. IBC I GO TO 120 
im=( jbc-ii*ms-jbc* uec-i )/2 ♦ IRC 

X 3= X I 4, IM) 

X4=X (3. I M I 
GO TO 121 

120 X 3= X ( 3, IM ) 

X4=X(4,IM> 

121 Xl=X!l,IMJ 
X2=X (2, IM) 

X5 = X I 5 , I M) 

COMPUTF THE UPPER DIAGONAL ELEMENTS OF THE MASS AND STIFFNESS 
MATRICIES 


CONTRIBUTIONS OF SHELL 

IF I J .LT. I ) GO TO 112 
T F ( J .LE. 0 .OP. I.LE. 0) GO TO 114 
f SUBMATRIX A 

AM I, J)=AM I.J >*Sl*R 1 1) *X1*I S2*R I 2) *S3*R I 3) ) 4ABN* X2 
C. SUeMATRIX N 

AMI I , Jl =AMt I » JI + PC *R ( 1 1 * X2 
C SUBMATRIX B 

114 AK(IN,JN)=AK<IN,JM*S1*ABN*RIB)*X5*IS2*R(<))*S6*RI2))*X2 
C SUBMATRIX 0 

AMI !N,JNJ*AMt IN, JN|*PC*R(9I*X5 
C SUBMATRIX C 

AKI INN, JNN I = AM INN,JNN)*t SI *R 18) ♦ S5*R 14 1 ) 6X5* S5* I RI l I *X 1 ♦ I ABN2-AN2 
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non 


1-BN2 »*RIA)*X5 )-S7*R(8I*I8N2*X3aaN2*XA )aS8*ABN*R( 2)*X? 

C SUBMATR l X S 

AMI INN, JNN)«AMI INN, JNN I APC *R I 1)*X5 

112 I F ( I .IE • 01 GO TC 113 

C SUBMATRIX 0 

AM !, JN) = AK< I,JN)aSA*BN*R( 5 )*X3-S2*AN**t fc)*X2 
C SUBMATRIX E 

AM I, JNN ) = AKI I, JNN UR I 51*1 S A* X3- S5* X 1 ) a S 3* ABN*R I ?) * X2 
C SUBMATRIX F 

113 AKI IN, JNN) =AM I N , JNN) ♦ AN*R ( 8 ) » I S 1 *X5 -S7* X A I »S6 * BN*R 1 2 » *X 2 
IF I IR ,E0. 0) GO TO 111 

I F I J .IT. M GO TC 115 
C SUBMATRIX B 

AM IN,JNI = AK( IN,JN)AS5*RV( 1I*X5 
C SUBMATRIX F 

115 AK( IN, J^N)=AK(lN,JN.'i|-S7*KVl 3 I* XA A $S*( oN*R V(l I ~ ( 1.0 -t)N2)*HVU) ) • 

1 X5 

C SUBMATRIX C 

AKI INN, JNN )= AKI INN, JNN) AS5*lAbX*nV( 11 ♦! ABNB- A N) • * Vl 2 I ♦ I Ad NA- BN) * K V 
1(A) )*X5-S7* I 6N*RV 13 )*XJ«AN*KV 13 >*XA ) 

CCNTR I BUT IONS OF STRINGER 

111 IFILL .EO. 0) GO TO 106 
I F ( J .IT. I) GO TC 116 
IFIJ .LE. 0 .OR. I .LF. Cl GO TO 117 
C SUEHATRIX A 

AKII ,JI = AK 1 1 , J) a $T 1 1 I *X 1 
C SUBMATRIX N 

AMI I , JI = AMI I, J UST (29 )*X2 
C SUBMATRIX 8 

117 AKI IN, JN ) = AK I IN.JNUSTI 2 I * X l* ST I 7 A) * X2 

C SUBMATRIX 0 

AMI IN,JN)=AM( IN,JNItST( 30* X2*ST131 )*X5 
C SUBMATR IX C 

AKI INN, JNN) = AKI INN, JNN) AST I3)*X1 AST (73)*X2 
c SUBMATRIX S 

AMI INN, JNN 1= AMI INN, JNN ) a$T( 3 3 ) * X 2*1 ST I 3 A ) ♦ S T( 35 > ) * X 5 
C SUBMATRIX R 

lit AMI IN, JNN) =AM( IN, JNN 1 ♦ S T I 32 ) * X5 
C SUBMATRIX F 

AKI IN, JNN) =AK(IN,JNN)*ST(75 )*X2 
IF! 171 .FQ. 0) GO TO 107 
IFIJ .LT. II GO TO 118 
C SUBMATRIX B 

AKI IN.JN l=AKI IN. JNUSTI 5)*X1 
C SUBMATRIX 0 

AMI IN.JN) =AM( IfJ.JN) *• S T I 371* X2*ST 138 >*X5 
C SUeMATRIX C 

AKI INN, JNN) = AKI I NN, JNN) ♦ 1ST 17 ) AST (8 I )*Xl 
C SUBMATRIX S 

AMI I NN, JNN I a AM I INN, JNN ) ♦ I ST ( A 1 I a$T ( A2 ) I *X 2a S T I A 3 1 * X 5 

118 IF 1 1 .LE. 0) GO TO 119 

C SUeMATRIX E 

AM t , JNN 1 = AK I t , JNN I AST ( A I *X 1 
C SUBMATRIX P 

AMI I, JNN )= AMI I, JNN) ASTI 36)* X2 
C SUBMATRIX F 

US AKI IN,JNNI=AKIIN,JNN|aST(6)*X1 
C SUBMATR IX R . 

AMI IN, JNN) = AM ( IN, JNN) AST (39 ) *X2 aST I AO ) *X5 

1C7 IF! 122 .EQ. 0) GO TO 108 
IFIJ .IT. I) GO TO 520 
C SUBMATRIX 0 

AMI IN, JN ) = AM( IN, JN)AST( A5)*X5 
C SUBMATRIX C 
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A K ( INN, JNN) =AK( l NN, J NN I ♦ S T I I C ) • XI 
C. SUBMAT R I X S 

AM( I NN, JNN) *AMII NN, JNNI *ST UT )*X5 *ST (48 » *X2 

520 IF! I .IE. 01 GO TO 521 

C SUPMATBIX F 

AK ( I , JNN) =AK( I (JNNI ♦ STI 9 ) *X 1 
C SUBMATR 1 X P 

AM( l, JNN NAN I 1, JNN I *ST 1 44 ) *X2 
f. SUBMATRIX P 

521 AMI IN.JNN )=AM( IN,JNNI*ST( 46)*X5 

108 IFII71 .FO. 0) GC TO 109 

IFtJ .IT. I) GO TO l 22 
C SUBMATRIX B 

AKIIN.JNI = AK(IN, JM*ST (13)*X1 
C SURMATRIX 0 

AMI IN,JN)=AMI IN, JNI*ST(5l)*X2*ST( 52I*X5 
C SUBMATRIX C 

AM I NN, JNN ) = AK 1 1 NN, JNN I ♦ I STI 1 6) ♦ STI 1 T I ) * X 1 
C SURMATRIX S 

AM( INN, JNN) = AMI I NN , JNN ) ♦ 1ST I57I*ST(58) )*X2 ♦ I ST I 59 I *ST (601 )*X5 

122 IFI I .IE. 0) GO TO 123 

C SUPMATRIX 0 

AM I , JN I =AK I I ,JN)*ST( 11 1* XI 
C SUBMATR IX NN 

AMI I ,JN) =AM(I , JN)*ST (49 )*X2 
C SUBMATR T X E 

AK I I , JNN )= AK I I, JNNMSTI 12)*X1 
C SURMATRIX P 

AMI I, JNN ) = AM I I, JNNI* STI 5C)*X2 
C. SUBMATRIX F 

123 AKI IN, JNN)=AKI I N , JNN > ♦ I S T 1 1 4 ) *S T 1 15 ) I*X1 

C SUBMATRIX P 

AMI IN, JNN )= AMI IN, JNN ) ♦( ST 1 5 3 ) *ST I 54 )I*X2*( ST ( 55 ) ♦ ST l 56) I * X5 

ICS IF! I V2 .EQ. 0) GO TO 110 
IFI J .IT . II GO TO 124 
C SUBMATRIX R 

AKI IN,JN)=AK(IN,JN)*ST(2C)*Xl 
f. SUBMATR IX 0 

AM! I N ,JN) = A M { I N , J N ) * S T <6 3 )*X2 *ST (64 I *X5 
C SUBMATRIX C 

AKI INN, JNN) = AKI INN, JNN » ♦STI 23)*X1 
C SUBMATRIX S 

AMI INN, JNN )= AMI I NN, JNN) *1 ST I 69) *ST< 70) ) * X2* ( STI 71 1 *ST 172 ) l *X5 

124 IFI! .IF. 01 GO TO 525 

C SUBMATRIX 0 

AKI I, JN ) = AK I I , JN >*S T I 1 8 ) *X 1 
C SUPMATRIX NN 

AMI I , JNI=AMI I ,JN)*STI61I • X2 
C SUPMATRIX F 

AKI I , JNN) =AMI , JNM*ST I I9I*X1 
C SUBMATRIX P 

AMI l,JNNI=AMl!,JNN)*STI62)*X2 
C SUBMATRIX F 

525 AKI IN, JNN) = AKIIN,JNN )*I STI 21 )♦ STI 22 I >*Xl 
C SUBMAT 9 IX R 

AMI IN, JNN) =AMIIN, JNNI* I ST (65 l*ST 166 ) l*X2 *IST 167 )*ST (68) l*X5 

110 IF 1 1 2 1 .EQ. 0 .OR. I VI .EQ. 01 GO TO 106 
IF I J .IT. II GO TC 526 
C SUBMA TR I X B 

AKI IN,JN)=AK I IN, JN ) *ST( 24 )*X 1 
C SUBMATRIX C 

AKI INN,JNN)=AK(INN, JNN)*(ST(27)*ST(28))RX1 


I 
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C SU8MATRIX F 

526 AX ( I N» JNN ) =AK( I N , JNN )♦ I ST 1 2 5 ) *S T I 26 ) l*Xl 

CONTRIBUTIONS OF RING 

106 IF ( KK .60. 0) GO TO 93 
00 166 KR=1,NK ’ 

>Xl=XXXll ,KR,IM) 

XX2*XXX( 2,KR, INI 
IF1 J .IT. !> GO TC 151 
I F ( J .16 . 0 .OR. I . 16. 0> GO TO 153 
SUBHATR IX A 

AM I , J ) = AM I «J)MCR(KR, 1 I ♦ ABN2 *R 1 (KR , 1 ) »CR I KR , 2 1 I *ABN *R 1 1 KR , 51 1* 
1XX1 

C SUeMATRIX N 

AMU ,j) = AM(l , J|+(CR(KR,22I*RI (KR , 6 1 *CR I K R, 2 3 l *R l IXR , 7 1 *ABN) *XX l 
C SUBMATRIX B 

1 53 AXI IN,JN)=AK(IN»JN)-»-(CR(KR,2l*RI (KR , 3 ) ♦CR (XR , 3 ) *R I ( XR , 1 1 1 *ABN*XX2 
C SUBMATRIX Q 

AMI IN,JN)=AM( IN, JNI+ICRIKR, 22I*RI I KR , 8) *CR I KR ,26) *RI I XR ,7>)*XX2 
SUBMATRIX C 

AKI INN, JNN)=AK( INN, JNNI ♦ (CR I KR , 3 1 * ABN2*R I IXR ,1 I *CR I KR ,2 I *Rl IXR, 3) ) 
l*XX2*ICR(XR, 1 l*RI(KR,3)*CR(KR,2l 1*ABN*R H XR ,7))*XX1 
SL'BMATRI X S 

AMI INN, JNN) =A Ml I NN, JNN) ♦ICRIXR, 26)4 CR IXR, 231 >*R 1 1 KR,6 ) *XX 1* I CR (KR , 
1 26 ) * ABN*R I (KR,71»CR(KR,221*RI(KR«6)I*XX2 

151 IF t I .LE. 0) GO TO 152 

C SUBHATR IX 6 

AKI I , JNN I* AX I I, JNN) ♦ICR I XR , 1 ) *AN2 *R 1 1 XR , 2 ) ♦CR I KR , 2l)*ABN*RI I XR ,6) > 
1*XX1 

C SUBMATRIX F 

152 AKI IN,JNN) = AKIIN»JNN)*ICRIKR,3)*A8NB*R!(KR, 1 ) *CR I KR , 2) *AN*R II KR , 31 
1 )*XX2 

C SUBHATR IX R 

AM( I N, JNN) = AM (IN, JNN ) * 2 .0 E0*CR I KR , 26 ) *BN *R l I KR , 7 ) *X X 2 
IFIEIRIKR) .60. 0.0 ) GO TO 156 

IF | J .LT. II GO TO 155 
C SUBMATRIX B 

AX I IN, JN ) = AKI !N,JN)MCR(KR,5)*RI IXR,9)»CRIXR,6)*RI (KR»10)*CR(KR,7) 
1*RI IKR.ll I +CR IKR,8 )*RII XR ,121 1*ABN*XX2 
C SUBMATRIX 0 

AMI IN, JN ) = AM( IN, JNIMCRIKR, 30)*RI (XR, 15) *CR I XR , 2 8) *R I I XR , 16 1 ♦ 
l CRI KR ,32 ) *R I (KR,17I«XR(KP,33)*RMKR, 18) )RXX2 
C SUBMATRIX C 

AXI INN, JNN ) = AXI INN, JNN) ♦< I CR I KR , 51 *R I I KR , 5) *CR I KR , 6 1 * RII KR , 1 0) ♦ 
1C.RI KR ,81 *RI IKR,12ll*ABN2»CR(KR,q|*RHXR, 1 1 I * I AN 2*BN 2 1 I *X X 2* (CR IXR , 
2 10 1 *R 1 1 XR , U*ABN2-CR(KR,6)*RI(XR,2)*IAN2»BN2)»(CRIKR,26>*PI IKR,5>- 
3CPIKR.27 )*Rl (XR ,6)1 *ABN )*XXl 
C SUBMATRIX S 

AMI INN, JNN J = AM( INN, JNN I ♦ I CR I XR , 30 ) *R I IXR ,6) ♦ CR I KR ,36 |* ABN*R I (KR ,7 I 
1 l*XXl + (CR(XR,30l*RI(KR, 15UCRIKR, 33 1 *R 1 1 KR , 1 8 ) ♦ CR I KR , 32 ) *R 1 1 XR , 1 7) 
2)*ABN*XX2 

155 IF I I .16 . 0 1 GO TO 156 

C SUBHATR IX 6 

AKI l , JNN) =AK( l , JNN) -ICR I KR ,6 > * ABN2* R I (KR,l I »CR I XR ,25 ) *ABN*R 1 IXR, 5) 
l )*XXl 

C SUBMATRIX P 

AM( I , JNN I *A Ml I ,JNNI-(CR(KR,281*RI (KR ,6 ) ♦CR I XR ,2 9 ) *ABN*R I (XR, 7 ) ) * 
1XX1 

C SUBMATRIX f - 

156 AKI IN, JNN)* AK( IN, JNN ) ♦( (CR(KR,8I»RI I XR , 1 2 1 *CR I KR ,5) *R I (KR.R) ♦CRIKR 
l ,6 I *R I IKR.10I I *ABNBM ABN8*AN ) *Cfi I XR, 9 1 *R I (KR, ID )*XX2 
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C SLBMATRI X R 

AM( IN, JNN )= AM ( IN, JNN I *2 .OEO*BN*( CR ( KR , 31 ) *R ! ( KR , 16) *CR( KR *30I*R1 1 K 
IR.15)*CR(KR,33I*RI (KR, 1 8 > *CR (KR ,32 ) *R H KR l 7 ) )*XX2 

154 IF(E2R(KRI .EO. 0.0 > GO TO 157 

!F(J .IT. II GO TO 158 
C SUBMATRIX B 

AK ( IN,JNI=AK( IN,JNI*(CR(KR,11»*RI(KR,2)*CR(KR,12)*RI ( KR, 19 ) *CR (KR , 
113I*RI IKR,20I )*ABN*XX2 
C SUBHA TR IX 0 

AM ( IN,JN)=AM( IN,JN>»(CR(KR,35I*R(6I+CR(KR,37)*R(7)*CR(KR,38I*R(2>) 
l*XX2 

C SUBMATRIX C 

AK ( INN, JNN) = AK I INN, JNN I ♦( ( C« (KR , 121 =»R HKR , 19 5 *CR S KR , 1 7> *R I ( KR ,201) 
1*ABN2*(CR(KR,16)*RI ( KR ,201 ♦ CR ( KR ,1* )*Rl (KR»2 S IM AN2*BN2 I l*XX2 
C SUBMATRIX S 

AMI INN,JNN)=AM(INN,JNN) *CR (K R , *0 I *R I (KR , A I *XX l ♦( CR ( KR , AO 1 *R ( 21 *CR ( 
1KR,37)*R(7) l*ABN*XX2 

158 IF( I .LE. 0 I GO TO 159 

C SUBMATRIX P 

AM( I , JNN)=AM( I , JNN) -CR ( KR , 3 5) * R I ( KR , A I*XX1 
C SUBMATRIX F 

159 AK ( I N , JNN) =AK(IN,JNN)*(CR(KR»14 ) *Rl (KR, 2 I ♦( ABN8* AN I ♦ ( CR (KR, 15)*RI( 
IKR»20)*CR(KR,12)*RI(KR,19)I*ABNB*CR(KR,16I*R! ( KR „ 20) * AN) • XX2 

C SUBMATRIX R 

AM{ IN, JNN) =AM(IN,jNNI*(CR(KR,35 l*R(6 1 *CR (KR, 39I*R(2I ♦2.060*CR(KR , 

1 37 ) *R ( 7 ) )*BN*XX2 

157 IF ( I R .EO. 0) GO TO 144 

IF ( J .LTi I I GO TO 160 

!F(J .LE. 0 .OR. I .IE. 0) GO TO 161 
C SUBMATRIX A- 

AK( I , J) =AK( I ,J)*CR(KR, 1)*(ABN*RI (KR ,2 1 ) ♦ A8NB*R I ( KR ,22 ) ♦ ABNA* R! (KR, 
123 I )*XX1 

C SUBMATRIX B 

161 AM IN , JN) =A M IN,JNI+CR( KR,3) *(RI ( KR , 2 1 ) ♦ EN*R I (KR ,22 I *AN*R I (KR, 23 ) I 
1 *XX2 

C SUBMATRIX C 

AK( INN, JNN) =A K( INN, JNN) ♦CR(KR,3)*( ABN*R l 5KR,21 I »ABN8*R I ( KR, 2 2 )♦ ABN 
1A*RI (KR,23I)*XX2 

16C IF ( I .LE. 0) GC TC 162 
f SUBMATRIX F 

AK ( I , JNN)= AK ( I , JNN ) ♦CR ( KR, 1 ) *AN *R I ( KR , 24 1 *XX 1 
C SUBMATRIX F 

162 AK ( I K, JNN )= AM IN, JNN ) *CR ( KR , 3 ) •( BN2*R 1 1 KR , 22) *Rl I KR ,23 > 

l * ABN* BN* R I (KR,2ll)*XX2 

IF( E IR( KR) .FQ. 0.0 I GO TO 163 

IF ( J .LT. I) GO TO 164 
C SUBMATRIX B 

AK ( IN, JNI=AK( IN, JN I ♦ ( CR ( KR , 5 ) •( RI (KR ,43>+BN*RI IKR.35) *AN*R! (KR.36) 

1 )*CR(KR,6)*(RI (KR,26 I *R I ( KR , 27 ) *BN*( R l IK R, 37)*R I (KR , 39) ) +AN *( R 1 ( KR 
2,3B)*RI ( KR,40) ) )♦ 

3 CR(KR,8)*(RI(KR,28)*RI(KR,29I«-AN*RI(KR,42)*BN*R|(KR, 

441 ) I *CR ( KR,9 )*( AN*RI (KR, 32 I *BN*R I (KR, 31 ) )>CR(KR, 1« ) *R I ( KR , 30) *CR ( K 
5R,19)*(AN*RI(KR, 34)*BN*RI (KR,33) ) )*XX2 
C. SUBMATRIX C 

AM I NN, JNN) =AK( I NN , JNN) ♦ ( CR (KR, 1 0 I* ( R I (K R, 2 1 l<=ABN*ABNB*R I (KR , 22) ♦ 

1 ARNA*R I ( KR , 23) l-CR(KR,4) *(BN*R( ( KR,25)*AN*RI (KR.24I ) )*XX1 ♦ ( (CPIKR, 
251 *R I (KR,43)*CR(KR,6 I ♦( R I ( KR , 26 I *R l ( KR , 2 7 1 ) *CR ( KR , 1 8) *R( (KR,30)*CR 
3(KR , ei*(RI( KR ,29| +R( (KR ,281 I I* ABN* C R (KR ,5 I * ( C0NB*R I (KR , 35 ) MBNA * 
4RI (KR, 36 ) )*CR(KR,6)*(ABNB*(RI(KR , 37 1 *R I ( KR , 39) l oABNA* ( RI ( KR ,38) ♦ 
5RIIKR.40) ) )*CR(KR, 8) *(ABNB*RI(KR ABNA*RI ( KR , 42 I I *CR » KR , 1 9 ) • ( AB 

6NB*RI(KR,33I*ABNA*RI(KR,34) I *CR (K R, 9 ) *( BN*R I (KR , 32 1 *AN*fi l ( KR ,31) ) ( 

7 * X X2 
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164 IP ( ! .IE. 0> GO TO 165 

C SUBHA TR IX E 

AK ( I, JNN )=AK ( I , JNN l-CR ( KR ,4)*(ABN*Rl ( KR , 21 I ♦ ABNA* Rf ( KR, 2 3 » ♦ ABNB* 
IRl (KR,22 ) )*XX1 
C SUBMATRIX P 

165 AK ( IN# JNN )= AK I IN* JNN ) ♦( (CR(KR,5)*RI(KR,35)* 

* CR(KR,6)*(RI(KR,37)*Rl(K 

1R«39))*CR(KR,8)*RI(KR»41>*CR(KR»19)*RI(KR»33) > *BN2 ♦ ( CR( KR,9 ) *Rl (KR 
2,32>*CR(KR,5)*RI (KR, 36> *CR(KR, 19)«« I ( KR , 34)*CR( KR ,6)*(Rl (KR,38)*Rl 
3 (KR,40) ) *CR ( KR, 8 ) *RI (KR.42 I I*ABN*(CR (KR, 5 »*R KKR ,43 •♦CR (KR, 61 •( R1 ( 
4KR, 26l*Rl(KR, 271 ) *CR ( KR , 8) • ( R H KR ,2 8 ) *RI ( KR.29I I *CR(K«, 18 )•«! IKR,3 
50) )*BN*CR(KR,9I*RI(KR,31 ) )*XX2 

U3 IP ( E2R ( KR) .EC. 0.0 » GO TO 1*4 

IF( J .LT. || GO TO 166 
C SUBMATRIX 8 

AK< IN, JN)=AM IN, JN) ♦ (CR ( KR , 1 2 ) * ( RIIKR,** )*RI (KR,*6 I ♦8N*R I (KR ,49 > 
l+AN*RI(KR,50)>»CR(KR,17)*(RI(KR,45)*BN*(RI(KR,47>«ftHKR,53))*AN*( 
3RKKR.4B >*RI (KR,5*I I I +CR (KR, 14). * ( BN*R KKR, 2* KAN *R I (KR , 25) > ♦CRl KR , 
4201 *(BN*RI ( KR »51)*AN*R1 ( KR.52) > >*XX2 
C SUBMAT R IX C 

AM I NN, JNN) = AM INN, JNN) ♦ (CR (KR, 12 )*( ABN*( R I (KR, 44 KR I (KR , 46 1 I ♦ABNB 
l*RI (KR, 49) *A8NA*RI(KR,50)KCR(KR , 17)*RI ( KR ,45)* ABN*C R (KR , 20 > • ( ABNB 
2*RI (KR.51 KABNA*Rt(KR,5 2 )KCR(KR, 16 ) ♦( ABNBM R I (KR , 5 3 1 *R ( ( KR » 47) I ♦ 
4ABNA*(RI(KR ,54KRt (KR.48I ) ♦ BN* ( R I ( K R ,48 I *R I (KR, 54 ) ) ♦ AN* (R KKR ,47 )♦• 
5RI (KR,53 I ) ) *CR( KR, 14) *( AN*R K KR , 24) *BN*RI (KR ,25) I KXX2 
C SUHMATRIX P 

166 AK ( IN, JNN) = AK( I N,JNNK( ( CR ( KR , 1 2 ) *R I { KR ,49 )*CR( KR, 20 ) *R I ( KR, 5 1 ) ♦ 
1CRIKR, 16)*(RI (KR,53)*RI(KR,4?)) ) *BN2*(CR ( KR , 17) M R I ( KR, 4BKRI (KR, 

2 541 ) *CR ( KR, 12 KRI (KR ,50 ) *CR ( KR , 1 4 KR I (KR , 25 ) *CR (K R, 20 KR I (KR , 52 ) ) 
3* ABN ♦ ( CR ( KR , 12>*(RI(KR,44KRI ( KR ,46) ) *CR ( KR , 1 7) • R I ( KR.45 ) KBN*CR ( 
4KR, 16 K ( RI (KR.47 )*RI(KR,53) J*CR(KR, 14KP KKR, 24) ) *XX2 

144 CONTINUE 
93 CONTINUE 
92 CONTINUE 
SI CONTINUE 
90 CONTINUE 

COMPUTE THE LOWER DIAGONAL ELEMENTS OF THE SYMMETRIC MASS AND 
STIFFNESS MATRICES 

C 

00 2 CC2 1=2, MN3 

11=1-1 

DO 2002 J = l ,!l 

AM t,JI=AK( J , I ) 

2002 AM( I, J ) = AM( J, I I 

IF (NWK-1 12010,2011,2012 

2C12 WRITE( 7,65) ((AM I , J) ,J=1,WN3) ,1=1 ,MN3) 

2011 WRITE (6 ,2015) 

2015 FORMAT! 1HI ,30X ,1 7H$T IFFNESS M ATR IX , / ,3 IX , 17 ( IH= ) , // ) 

00 2013 1=1, MN 3 

WRITE (6, 201 4) (AK( |,J ),J=l,MN3) 

2014 FORMAT (5X,E15.8,1X,E15.8,1X,E15.8,IX,E15.8I 

2013 CONTINUE 

2010 IF(NWM-1 >2016,2017,2018 

2016 WRI TE( 7,65) ( (AM( I , J) ,J=1 ,MN3) ,1 =1 ,MN3) 

2017 WR I T E ( 6 , 20 19 ) 

2019 FORMATllHl ,32X,12HMASS MAT R I X, t , 33 X , 12 ( 1 F= ) ,// ) 

00 2020 1=1, MN3 

WRITE (6, 20 14) I AMI I,J ),J= 1.MN3) 

2020 CONTINUE 
C 

C EVALUATE THE EIGENVALUES AND EIGENVECTORS 
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non 


c 

2016 GALL EIGENO (AM ,AK ,MN3, KRRR , VECR ,LC ,*XXX,V„**C ,Z, E VR»F Vl.INDI C » 

CONVERT THE EIGENVALUES INTO FREQUENCIES IN HERTZ 

CC 2006 1*1, MN3 - 

IF (EV.RIT I. .LE. 0.0 1 GO TO 2006 

EVRI l ) = SORTIEVRM II/PI 2 
2006 CONTINUE 

IFINWEV-1I2021, 2022,2023 
20:23 00 2026 J= l,MN,3 

WRITE IT, 65 I LVECRT I.,J.», !*l,MN3:l 
2.0.26 CONTINUE » 

20 2:2 WRITE('6„2025) 

2025 FORM AT < l HI , 32X, 1 2HE IGENVECTORS, A, 33X, 12( IH» » , /;/»; 

TOT 2026 J=l ,MN3 . 

WRITEI6,2016MVECR( I,J», 1*1, MN3) 

2026 CONTINUE 

2C.2 V WRI TE! 6, .20271 

2027’ FORMAT ( l HI ,27X , 22HE I GENVALUE S IN HER TZ , A, . 2.0 X, 221 IH*| ,//» 

WRITE I 6, 201 6) IEVRIII ,I = 1,MN3). 

GO to 10000 
R9 WRITE <6, 100 1, 

liCO FORRATIr//, l.X,62H***R* ERROR ****•• MM I N = 0 I S POSSIBLEONLYWT 
1TH FREE-FREE, /, 32X, l 8HB0UNDART' CONOITION) 

CC TC 10000 
100C5 CONTINUE 
CALL EXIT 
ENO. 



U 
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*********************************************************** ********* 

* ’.'■'♦ 

♦ SUBROUTINE INTGRL ♦ 

♦ ♦ 

♦ PURPOSE , , ♦ 

♦ TO EVALUATE THE LONGITUDINAL INTEGRALS XI (II TO XtI5) WHEN ♦ 

♦ N , M0, AA, ANC NBC ARE FURNISHED TWOUC-H THE COMMON STATEMENT ♦ 

♦ * 


♦ USAGE 

♦ CALL INTGRL 

♦ 

♦ DESCRIPTION OF THE PARAMETERS 

♦ PM 

♦ MB MB 

♦ . AA LENGTH OF T PE SHELL 

♦ NBC COOE NUMBER OF THE BOUNDARY CONDITION WHICH IS 

♦ UNDER CONSIDERATION 

♦ PI 3.1*15926535 

♦ XIII) LONGI TUCINAl INTEGRAL- 


♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 


♦ SUBROUTINES required ♦ 

♦ NONE ♦ 

♦ ♦ 

♦ METHCC ♦ 

♦ the CLOSE FORM EXPRESSIONS FOR THE INTEGRALS WERE CBT A I NEC ♦ 

♦ FROM ROBERT P. FELGAP, JR., FORMULAS FOR INTEGRALS ♦ 

♦ CONTAINING CH AR ACT ER I ST I CFUNCT IONS OF A V IBRAT ING BEAM, ♦ 

♦ CIRCULAR NUMBER 1*. THE UNIVERSITY OF TEXAS, IPSO ♦ 

♦ ♦ 


SUBROUTINE INTGRL 
C I ME NS I ON A(5I« PI5I 

COMMON 0K(91 ,RI91 , ORV 15 I ,k v IS I ,RI <u » , KR1 Id ) , *2 UO ) , RR 2 HO ) , R iU I , 
1RR 3 1 2 ) , R At 5)»RR*{51.R5< 181, R RSI 181.R6I11I , RR6I l 1 1 ,PI , XK,AA,Xt<5) , 
2XRI2) ,E1 ,E2 ,N,NB t NBC,M, MB.NSA 
DO 8 1=1,5 
XIII 1=0 .0 
8 CONTINUE 

GO TO I 100, 200, 3CC, *C01 ,NBt 
CLAMPEO-FPFE 

ALPHAS AND BETAS CBTAINEC FROM FELGAR ANO YOUNG 
100 A I 11 = 0. 7340955 
A(2 1=1 .018*66** 

AI3I=C. 99922*50 
A{*)=1. 000033553 
A 15 1=0.9999985501 
B( 11 =1. 8751C*1E0/AA 
RI2 l = *.69409U3E0/AA 
et3 1 =7.R5*757*3E0/AA 
BI*1=1C. 99 55*07 AEC/AA 
PIS »=l*.13716839E0/A A 
IF I M - 51 15,15,25 
15 AM= AIM) 

PM=e IM) 

GO TC 30 

25 BM= ( 2. OEO* FLOATIMJ-I.OECIPPI /I 2 • OE 0* AA 1 
AM=1.0 
2C AB= AM*BM 

P*= BM*BM*BM *BM 
IF IMP - 51*0, *0,50 
AC AMB = At MB I 
PMB=B(MB » 
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GO TO 20 


«C AM8 = 1.0 

8Mfl=(2.0E0* FL0AT(M8I-1.0E0I*PI/I 2.0EO*AA» 

20 A8R=AMB*BMB 

P4B= bmb*bhb*bmb*bmb 
IFIM-MB) 60,65,60 
f. M = MB 

65 XIII J=84*AA 

XI (2 I =AB* 12 .0 ♦ AB*AA I 

XI l 3 I =AB*I 2.0 -AB*AA) 

XI(4I=X II 31 
X 1 1 5 I =AA 
. , BE TURN 

C M NOT = MB ' 

60 BBM=BM*BM8 
MPMB = N«-MB 
A3B» AMB *BM*8M*8M 
AB3 =AM*RM8*8MB*BMB 

XII 2I=4.0E0*B8M*< (- I .OE 0) **MPM8* ( A3B-AB3 1 -88M* I AB-ABBI 1/ ( B4-B48I 

X 1 1 3 1=4 .OEO*BM*BM* I ABB- AB ) *1 I I - 1 .OE 0 I •*MPMB » *BM*RM*BMB*BMB> / 
1IS4B-B4) 

XII 4 1 = 4 .OEO*BMB *8M8 *1 AB-ABB ) •( ( I - 1. OE 0) * • MPMBI * 8MB* BM6»8M*BM I / 
l (B4-B4B) 

RETURN 

C FREELY SUPPORTED 

2C0 IF|M-M8)70,T5,70 
C. M b MB 
75 PM=M 

Xl(2l=BM*BM*PI*Pt/AA 

XI I 31 =- XI (21 
XI(4I=XI(3I 

XI m=xn2)*xl(2 l/AA 

XII 5)=AA 
7t) PFTUPN 

r clampfo clamped 

3C0 A I 1 ) =0. 08250221 58 
A (2 J=l .000777311 
A( 31 =0. <JS<J<J66450I 
A( 4 ) = l . COOOO 1450 
A (5 I =0.9999999373 
B( II =4. 730O40BF0/AA 
6(2 l=7.BS32046E0/AA 
PI3 I =10 .9956078E0/AA 
B( 41=14. 1371655F0/AA 
PI 5 1 = 1 7 . 2787596F O/AA 
IFIM - 511 ,1,2 
I AM=A (MI 
PM= B I M I 


GO TC 3 

2 BM*( 2. OEO* FLOATIMI*l.OEOI*PI /I2.0E0*AAI 
AM= 1 .0 

3 IF! MB - 514,4,5 

4 AM8= A I MB I 
BM8= P (MB I 
GO TO 6 

5 AMB=1 .0 

PMB= (2 .OEO* FLOATIMBI ♦ l .0 EO I *P 1/ 12 .OEO*A A I 

6 AB=AM*BM 
ABB= AMB*BMB 


84= BM*BM*BM*BM 

B4B=BMB*BMB*BMB*BMB 

BB2=BM*BMB*BM*8MB 
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lF<w-M8)80,B5,80 
C M = MB 

85 XK1) = 84*AA 

XI <2)=AB*(AB*AA-2.0 > 

XI ( 5 ) = A A 
GC TC 90 

C M NOT « MB 

80 X 1(2 )=4.0FO*BB2*(AR8-AB>*<(-i.OEO)**(M»MB)*l.OEO) /(BA -BAB) 
9C XII3I-XII2I 
XI(4)=Xl(3) 

RETURN 

EPEE-FREE 

ALPHAS ANO BETAS OBTAINED FROM FEIGAR AND YOUNG 
A 00 All )=0 .9825022158 
A ( 2 1 =1. CCC777311 
A ( 3 ) = 0. 999966450 1 
A (4 ) = 1 .00000 1 A 50 
A( 5) = C. 9999999373 
Bf I 1 = 4. 73004 C 8E 0/A A 
H ( 2 1=7. 8532046 EO/AA 
B< 3) =1 C.9956078F0/AA 
P ( 4 )=14.1371655 C 0/AA 
P ( 5 I =17 ,27875 00 EO/AA 
«M1=M-1 
V 9M I = MB- 1 
P2=N«-MB-2 

I E ( M .LT. 2) GO TO 1730 
IF(P - 61 715,715,725 
715 AN = A(MMU 
BM=H(MMl) 

GC TC 730 

725 B M = ( 2 • OE 0* ELCA T ( MMl I *1 .0 EO ) *P\/ 1 2 .0 EO* A A ) 

AM= 1 . 0 

730 A B=AP*BP 

P4=BM*BH*BM*BM 

1730 IFIMB .IT. 2) GO TO 1720 
I F ( M B - 6)740,740,750 
74C AMB=A(MBM1) 
e M B=B(MBHl ) 

GC TC 720 
75C AMR = 1.0 

PMR= ( 2 .OEO* FLOAT ( MBM 1 I ♦ 1 . OE 0) *P I /( 2. CF 0* AA ) 

72C A)P=AWB*B*B 

P4B=BMB*BMB*BMB*RMB 


1720 CCNTINUE 

Sl=(l.C - ( - l . 0 l**YBYl) 
S 2= ( 1 .0 ♦<-1.0 

53 = ( l .0 -1-1 .0 >**MMl ) 

54 = ( 1.0 ♦(- 1.0 )»*M"1 ) 

S5=(l.O ♦<-1.0 1 * *M 2 > 

IF(MNl) 150,250,350 

150 IF(PB .NF. H) GO TO 151 


C P=0 M8=0 

XI 15) =A A 
RETURN 

151 IF(Me .LT. 2) RETURN 
C M = 0 MB GRFATFP THAN OR =2 

X 1(4 ) = 2 .0 EO * ABB 1 * S 1 
OFTURN 

25C IF ( **B . . NE . M) GO TO 251 
C P=1 MR = 1 


XI (2) =1 . OEO/ A A 
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XK5I-AA/12.0 

RETURN 

251 IE( **8 .IT. 21 HE TURN 
C M=1 HR GREATER THAN OR * 2 

Xl(2)=-2.0E0*S2/AA 
XI ( 41=T 2. OF O/AA-ABBl *S2 
RETURN 

350 IF ( BBN1 ) 450 *550 #650 
G M GREATER THAN OR *2 MB*0 

45C xil31=2.OE0*AP*S3 
RETURN 

C M GREATER THAN OR =2 MB* l 

550 XII21=-2.0E0*S4/AA 

Xt I 3)=S4*< 2.0EO/AA-AB1 
RETURN 

C “ AND MB GREATER THAN OR *2 

650 IFIM ,NE. MB) GO TO 651 

C M = MB 

XI(ll*AA*B4 

XI <2I=AB*IAB*AA46.0EC) 

XI (3 )=AB*<2 ,OEO-Ae*AA) 

XI I 41 *XH 3) 

X I ( 5 )= A A 
RETURN 

G M NOT = MB 

651 AB3= AM*RMB*BMB*BMR 
A3B= AMB*BM*BM«-BM 

XI ( 2 1=4. CEO*BM*BMB* { AB3-A3B) *55/ IB4.B-B4 I 

X 1 1 3 1=4 .0E0*B4*C ARB-ABI *55/ 1 B4B-B4) 

XI (4) =4.0F0*B4B* IAB-ABB14S5/ IB4-B4BI 
RETURN 

END 
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SUBROUTINE XX 
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r 
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♦ 

PURPOSE 


4 

c 

4 

CALCULATE 

THE VALUES OF XRtli ANO, XRI2) AT A GIVEN VALUE OF 

4 

f 

4 

X * XK 


4 

r 
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4 

r 

4 

USAGE 


4 

r 

4 

CALL XX 


4 
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4 


4 

DESCR IPT ION 

OF THE PARAMETERS 

4 

r 

4 

o I 

3. 1415926535 

4 


4 

*K 

THE VALUE OF X AT MHICH THE FUNCTIONS XRI1) AND 

4 


4 


X R I 2 1 ARE TO RE EVALUATEO 

4 

r 

4 

A A 

LENGTH CF THE SHELL 
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4 

NBC 

THE CODE NUMBER OF THE BOUNOARY CONDITION UNDER 

4 

r 
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CONS IOERAT ION 

4 

c 

4 

M 

M 

4 

l 
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MB 

MB 
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c 
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r 
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REHARKS 
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THE INPUT 
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r 
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COMMON STATEMENT 
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SUBROUTINE XX 

C THIS SUBROUTINE EVALUATES THE VALUES OF Kill, »»I2I AT A GIVEN 

C VALUE OF X 

DIMENSION A ( 51 ,B( 51 

COMMON OR (9 ) » K (9 ) ,ukV l3),KV(5),Kl(o)»kkl(<il,K2I 10) ,Kk <:( 10) ,K i( c) , 
1RR 3(21 » R5 ( 5 I « RR5 1 5 ) , R5 ( 1 8 ) , RR5 ( 1 8 )',R6 I 1 1 ),RR6Ill I , P I , XK , AA, X I ( 5 ) , 
2XR ( 2 I » E 1 , F 2, N,NB,NBC,M,MB,NSA 
XPtl)=0.0 • 

XRI2)=0.0 

CO TO < 100, 200, 300, 500) , NBC 
C CLAMPFO-FREE 

C ALPHAS ANO BETAS OBTAINEO FROM FELGAR ANC VCUNG 

100 A(1)=0 .7350955 
A( 2) =1. Cl 856655 
AI 3 ) = 0. 99922550 
AI4|=l. 000033553 
AI 5) =0. <3999985501 
R( 1>=1.8751051EC/AA 
P(2 )=5.69509113E0/AA 
B<3I = 7. P5575753E0/AA 
e< 4 ) = 10.99554C74EC/AA 
P 1 5 ) = 15. 13716839 EO/AA 
IF I M - 51 15,15,25 
15 AM= A I M ) 

RM=R<M) 

GO TO 30 

25 ek=I2.0E0* FLOATCM)-1.0EO)*P!/I2.0EO*AA) 

AM= 1 . 0 

3C I F ( MB - 5 ) 40, 4C, 5C 


I 
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40 AMB=A< MB) 

BMR=B( MB ) 

GO TO 20 
•5C AMB=1.0 

BM8 = < 2 . 0£ O* FLOAT(*6)-1.0E0)*PI /I2.0E0*AA) 

20 eMX=8M*XK 
BMBX =8MB*X X 
EB = E XP( BMX ) 

EBB= EXP (BMBX) 

?* SI N( BMX) 

SB = SIN(BMBX) 

C = CCS (eMX ) 

CB= COS(BMBX) 

SH=(EB-1 .OFO/EB I/2.0E0 
SMB = ( EBB-l .OEO/ EBB ) / 2 .0 EO 
CH=( EB*l.0E0/EBI /2.CE0 
CHB = ( EBB *1 .OEO/ EBB)/ 2. OEO 

XR II) =BM*B**B* (SM*S-AM*(CM-C))*ISHB»SB-AMB*( CMB-CB) ) 
XR( 2 ) * ( CH-C - AM* ( SH-S) )*( CHB-CB-AMB* ( SHR-SB) ) 

RETURN 

C FREELY SUPPORTED 

2 CO BM=M 

RMB=MR 
PICA=PI/AA 
XM=P I OA *BM*XX 
XMB=PI0A*8MR*XK 

XR( l)=2.CE0*BM*BMB*PI0A*PI0A* COS(X“)* CCS(XMB) 

XR ( 2 )= 2 .OEO * S IN ( XM ) * SIN(X»(B) 

RETURN 

C CLAMPED CLAMPED 

<00 A( l NO .9825022158 
A ( 2 ) =1 .0 00777311 
A( 3)=0. 9999664501 
A ( A ) = 1 .00000 1450 
A ( 5 ( =0.9999999373 
B( 1) =4. 73OO4C0FC/AA 
e( 2 ) = 7.853204f>E'VAA 
8(3) =10. 99 5 60 7 8 e 0/ A A 
B( 4)= 14. 1371655F0/AA 
P(5)=17.278759*E0/AA 
IF(M - 5) 1,1 ,2 

1 AM= A ( M ) 

PM=B(M) 

GC TO 3 

2 BM= ( 2. OEO* FLOAT(M)*l.0EC)*PI/I 2.0EO*AA) 

A*=l .0 

3 f F ( M8 - 514,4,9 

4 AM8= A ( MB ) 

BMB= B (Mfl ) 

GC TC 20 

5 AMB= 1.0 

eM8=(2.0EO* FLOAT (MB) »1 .OEO ) *P I / ( 2 .OEG*AA ) 

GO TO 20 
C FREE-FPEE 

400 A ( 1 ) =0 . 9325022 1 58 
A(2)=l. 000777311 
A ( 3 NO. 9999664501 
A (4 ) = 1 .000001450 
A{ 5) =C. 9599999373 
e( l 1=4 . 7300 408E0 /A A 
8(2) =7.8532046F0/AA 
B( 3 > = 10 • 99 56C 7 8E C /A A 
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P<4 I = 14.1371M5E0/AA 
B( 51=1 7.27875965 O/AA 
MM 1 = M- 1 
PBP1=M8-1 

IF C P* .LT. 2) GO TO 1730 
IFIM - 6) 715,715,725 
715 AP=A(MM1> 

PM=fl(HMl) 

GO TO 1750 

725 PM= ( 2 .0 EO* FLOAT (PM l l»1.0E0l*PI/(2.0E0*AA) 
AM= l . 0 

1740 BPX=BP*XK 

S= SIMBMXI 
C= COS(flHX) 

EB= EXP(RMX| 

SH = (FB-1.0E0/EB1 /2.0E0 
CH= ( E B ♦ 1 .OEO/EB 1/2. OEO 
1730 IF(KP .LT. 2 I GO TC 1720 
l F ( M B - 61 740,750,750 
740 AMB= A ( MBH1 1 
BMR=>B< XB^l 1 
GO TO 1760 
750 AMB= 1 .0 

B M B = ( 2 • OEO* FLOAT (PBMl 1 *1 .OEO 1*P !/ (2 .OEO*AA I 
1 76C RMBX=BMB*XK 

SH= SIMRMBXl 
CB= C0S1BMBX1 
FBB= EXPtBMBXl 
SHB=(EBP-1 .0 EO/ EBB 1/2 .OEO 
CHR=<FBB*1.0EO/EBBI/2.0EO 
1 7 2C IF! HM1 » 150,250,350 

150 IF ( P8 . NE. PI GC TO 151 

0. M = 0 MR =0 

XR ( 2 1= 1 .0 
RETURN 

151 IF 1 MR Mil 152,152,153 

C p=0 M P= 1 

152 XRI21 =XK/AA-0.5 
RETURN 

C P=0 MR GREAT FR THAN OR -Z 

153 XR (2) =CHB*CB-AMB* (SHB*S PI 
PFTURN 

250 IF (PBMl ) 251 ,252,253 

C M=l MB=0 

251 XR ( 2 1= XK / A A - 0 . 5 
RETURN 

C M=1 MB = l 

252 XR< l 1 = 1 .OEO/ ( AA*AA> 

XRI2 1 =XK*XK/ ( AA* AAI *0.25 -XK/AA 
RETURN 

C P=l M 0 GREATER THAN OR =2 

25? XRUI =BP8*( SHB-SB-AMB*! CHB^Cei 1/AA 

XR( 21= (XK/AA- 0.5 1 *( CHB*CB-AMB* ( SHB^Sfil I 

RETURN 

35C IFIMBMl! 351,352,353 
C M GREATER THAN OR EQUAL TO 2 MB=0 

351 XR(2 >=CH*C-AM*(SH »S 1 
RE TURN 

C M GREATER than OR EQUAL TO 2 MB = 1 

352 XR(l >=8P*(SM-S-AP*(CHHCn/AA 

XR( 2I = ( XK/AA-O. 5 )*(CH*C-AM*( SH»S1 1 

RETURN 


37 



M ANO MR GUPATEB THAU OR «2 
^ 5 ^ XRI l )a8M*S«8*(SH-S-AN*(CH»C II •( SM8-SB-ANB* > 
XR<2 l*(CH*C-A*«(SH«SI »*( CHB+Cft- AM8*<St‘8»$Bt I 
RP TORN 
END 
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♦ ♦ 

♦ SUBROUTINE SHELHIX) ♦ 

♦ ♦ 

♦ PURPOSE ♦ 

♦ CALCULATE AT A GIVEN VALUE OE X THE VALUES OF THE INTEGRANDS ♦ 

♦ CF THE FIRST SET CF C IRCUMFE® ENT I AL INTEGRALS OF THF SHEU. ♦ 

♦ ♦ 

♦ USAGE ♦ 

♦ IT IS USED AS AN ARGUMENT OF THE NUMERICAL INTEGRATION ♦ 

♦ SUBROUTINE ‘GAUSS* ♦ 

♦ ♦ 

*■ INSCRIPTION OF THE PARAMETERS ♦ 

♦ X THE INPUT VALUE CF X AT WHICH THE INTEGRANDS HAVE *■ 

♦ TO RE EVALUATED *• 

♦ OR (II INTEGRANOS ♦ 

♦ N N ♦ 

♦ M3 NR ♦ 

♦ NSA 0 WHEN COMPUTING THE SYMMETRIC MODE SHAPES ♦ 

♦ 1 WHFN COMPUTING THE ANTI SYMMETRIC MODE SHAPES ♦ 

♦ ♦ 

♦ QEMARKS ♦ 

♦ EXCEPT FOR X ALL THE INPUT AND OUTPUT PARAMETERS ARE ♦ 

» COMMUNICATED THROUGH THF COMMON STATEMENT ♦ 

♦ * 

* FUNCTION SUBROUTINES REQUIRED ♦ 

♦ RSHL ♦ 

« ♦ 


SUBROUTINE SHELL1 IX) 

CCMMCN 0RI9»,RI9», URVI5 l,RV( SI.KUtt >,RRH8),*2I 101 ,RR2( 101 ,R 3< 2) , 
10R3I2I,RA(S),RRAI5I,R5I 18I.RP5I 18I.R6I lD.RRtl 11) ,P! , XK,AA,XI( 5) , 
?«R<2» ,E1 .F2 ,N f N8,NBC,M t Me,NSA 
«N=X*N 
X NR= X*NB 
SR =R SHL ( X) 

1FINSAM, 1,2 

1 f.N* COS< XN» 

SN= SINIXN) 

CNfl = COSIXNB) 

SNR * SIMXNRI 
GO TO 3 

2 CN= S1NIXNI 
SN s CCSIXNI 
r.NR= SINIXNBI 
SNR= COSIXNB) 

1 DR ( 5 I =CN*C NB 
OR I 6 »=SN*SNB 
CR(1 ) = S P*CR 15 I 
DR ( S ) = SR *PR I 6 > 

CR 12 )=0R 16 >/SR 
DRIB) =0R 16 ) /SR 
DR ( 7| = 0R< 21 /SR 
DR (? I =0R ( Y | / SR 
DR I A ) =0R < B) /(SR*SR) 

RE TURN 
FNC 
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♦ SUBROUTINE SHELL2IXI ♦ 

♦ ♦ 

♦ PURPOSE ♦ 

♦ CALCULATE AT A GIVEN VALUE OF X THE VALUES OF THE INTEGRANDS ♦ 

♦ OF THE SECOND SET OF Cl RCUMFERENTIAL INTEGRALS CF THE SHELL ♦ 

♦ 1 ♦ 

♦ USAGE ♦ 

♦ IT IS USED AS AN ARGUMENT CF THE NUMERICAL INTEGRATION «• 


▲ 

SUBROUTINE 'GAUSS* 
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DESCRIPTION 

OF THE PARAMETERS 
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4v 

X 

THE INPUT VALUEOF X 

AT WHICH THE 

integrands have 
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♦. 


TO BE EVALUATED 




4 

♦- 

DRV ( I I 

INTEGRA NOS 



1 

♦ 
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N 

N 
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NB 

NB 
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NSA 

0 WHEN COMPUTING 

THE 

SYMMETRIC MODE SHAPES 

♦ 
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l WHEN COMPUTING 
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ANT, ISYMMETR IC 

MODE SHAPES 
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REMARKS 
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♦ 

EXCEPT FOP X ALL THE INPUT 

AND 

OUTPUT PARAMETERS ARE 

4 


COMMUNICATED THROUGH THE COMMON STATEMENT 


♦ 

♦ 

♦ FUNCTION SUBROUTINES REQUIREC 

> R SHI 

♦ RRRT 


♦ 

♦ 

♦ 

*■ 

*■ 

♦ 


SUBROUTINE SHFLL2IX 1 

CCMMCN DR<9),R(9),URV(S.),kV(SI,R1(B),KkL(8) ,R 21 10) , K« 2U 0) . R i U ) . 

1RR3 (2),R4( 5 I.RRAI 5 1 ,R 51 1 81 , RP 5< 1 81 , R 61 1 1 1 ,RR6(ll ) ,P I » XK , A A , X I 1 5 I . 

2XPI2 I «E1 »E2 *N»N8*NBC«M»MB»NSA 
XN=X*N 
XNB=X*NB 
SR =RSHL ( X ) 

RT = RRR T ( X > 

PT2 - RT *RT 

SR 2=1. OEO /(SR* SRI 

IFINSA) 1, 1, 2 

1 SN= S I N I X N I 
SNB = SIMXNBI 
CNB= COSIXNB) 

CN= CCS I XN I 
GO TO 3 

2 SN= COS I XN I 
SNB= COSIXNB) 

CNB= SIN(XNB) 

CN= S IN ( XN I 
RT = — RT 

2 DRV! 1)=RT2*SN*SN8/SR 
UR VI 3) =R T*CN* SNB 
CKV(4)=0RV( 5)*SR2 
RETURN 
ENO 
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♦ ♦ 

♦ SUBROUTINE RINGHX1 ♦ 

♦ ♦ 

♦ PURPOSE ♦ 

♦ CALCULATE AT A GIVEN VALUE OF X THE VALUES OF THE INTEGRANDS ♦ 

♦ <-0F THE FIRST SET OF CIRCUMFERENTIAL INTEGRALS OF THE RING ♦ 
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♦ 
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♦ 

4 

4 

4 

4 
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USAGE 

IT IS USED AS AN ARGUMENT CF THE NUMERICAL INTEGRATION 
SUBROUTINE 'GAUSS* 

DESCRIPTION OF THE PARAMETERS 

El 7-OI STANCE OF THE SHEAR CENTER CF THE RING FROM THE 

MIOOLE SURFACE OF the SFFLL 

E 2 Z-OI STANCE OF THE CENTPCIO CF THE RING FROM ITS 

SHEAR CENTER 

X THE INPUT VALUE OF X AT WHICH TS-E INTEGRANDS HAVE 

TO BE EVALUATED 
Rl(I) INTEGRANOS 

N N 

NB NB 

NS A 0 WHEN COMPUTING THE SYMMETRIC MODE SHAPES 

1 WHEN COMPUTING THE ANT I SY MMET R IC HOCE SHAPES 
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4 

4 


4 REMARKS 4 


♦ EXCEPT FOR X ALL THE INPUT AND OUTPUT PARAMETERS ARE ♦ 

4 COMMUNICATED THROUGH THE COMMON STATEMENT 4 

4 4 

4 FUNCTION SUBROUTINES REOUIREO ♦ 

4 RSHl 4 

4 4 


SUBROUTINE RING! (XI 

COMMON OK (9 ),K< 9 J , JKV I S I. Kv I 51 , (U l ol .HR i( d) cl iOl , (W I lO J .Kdltt . 
1RRT(2»,RAIS»,RR«I 5I.R5I 18» ,RP5I 181 ,R6I ill ,RR6I11J ,PI ,XK,AA,XI(5I , 
2XRI2I.E1 ,E2 ,N,NB,NBC,M, MB.NSA 
XN= X*N 
XNB=X*NB 
SR=RSHL (X) 

RC= SR*-E 14F2 
I F ( NS A 1 1 1 li 2 

1 C N* CCS(XN) 

SN= SIN(XN) 

CNB= COSIXNei 
SNB= SINIXNBI 
CO TO 3 

2 CN= SIMXNI 
SN= COSI XN) 

CNB= SINIXNBI 
SNB = COS(XNB) 

3 CC=CN*CNB 
SS=SN*SNB 
RU4| =RC*CC 
Rl(8MRC*SS 
R1 (3 MCC/RC 
P1I 71 =SS/RC 

R1 1 2 ) = R 113) /RC 

Rl (6)*R1 (71/RC 

Pl( 1>*P l< 21 /RC 

Pll5)=Rl(6)/RC 

RETURN 

ENC 
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REMARKS 
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SUBROUTINE RING2IXI 




COMMON Uiil 9) i 

,RI 9) ,L)KV( 5) ,*v( b) ,Ki (a) ,Rk 1 (61 ,K2 1 iO J .iiK2 liO ! ,rO It 1 , 



IRR3I 21 ,R*( 5) 

,RRM 51 t R5( 1 8) , RR5I18I ,R6ll 1 ) ,RR6(lt > , <M , XK , A A, X ! (5 » , 



?XR(2 J, El, E2,N,NB,NBC,M,MB,NSA 
XK=X*N 
XNB=X*NB 
SR=RSHL (X » 

RC=SR*6l*e2 
! P( MSA ) 1, 1, 2 

1 CM COS (XN ) 

SN* SI M XNl 
CNR= COSIXNBI 
SN8= SINIXMRI 
GC TC 3 

2 CN= SINC XN I 
SM COS (XN I 
CNR = SI M XN8) 

SNR= COSIXNB) 

3 CC=CN*CNB 
SS= SN* SNB 
RC2=RC*RC 
RCSR »RC*SR 

R2I 10)=SS/RCSR 
P219 1= R2 < 10 i/SR 
P2(8> =RC*SS/SR 
R 21 7 ) = R 21 8) /SR 
R2 (6 )=SS/RC2 
R2I 51 =R2 (6) /RC 
R2( 3 1= CC/RC SR 


42 



R2 ( 1 I =R2 13 » /SR 
R ?( 4 1 *R 21 11 /RC2 
P2 ( 2 > = R2 <4 l/SR 
RETURN 
PNO 
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SUBROUTINE 'GAUSS* 




+ 

♦ 

description 

OF THE PARAMETERS 




♦ 

♦ 

El 

Z-0 1ST ANCE OF THE 

SHEAR CENTER OF THE RING FROM THE 

♦ 

♦ 


MIODLE SURFACE CF 

THE SHELL 


♦ 

♦ 

E2 

Z-DISTANCE OF THE 

CENTROID OF THE RING FROM 

I TS 

♦ 

♦ 


SHEAR CENTER 




♦ 

♦. 
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REMARKS 
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♦ 
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OUTPUT PARAMETERS ARE 
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COMMUNICATED THROUGH THE COMMON STATEMENT 
SUBROUTINES REQUIRED 


♦ FUNCTION 

♦ RSHL 

♦ 

SUBROUTINE RINGS IX> 

lUMMUN UR I VI ,K(9) ,ciRVI 5) ,KV (5) ,kl Id) , RKI (Si, hi (ID ) ,RR2< 10 >,K il 2 I, 
I PR3 12 I » RA ( S I » RR4 (5 I » RS I 18 I * RR 51 18 1 1 R6 1 1 1 I *RR 6( 1 1 1 *P I # XK , AA » XI I SI > 
2XRI2I ,E1,E2.N,NB,NBC,M,MB.NSA 
XN= X*N 
XNB=X*N8 
SR=R SHL C X) 

RC=SR*E1*E2 
I F ( NS A 1 1 «1 ,2 

1 CN= COSIXNI 
CNB= CCSIXNBI 
GO TC 3 

2 CN= SIN( XNI 
CNB= SINIXNB1 

3 CC =CN*CKB 

R3I 2 1 =CC /I RC*RC* SR I 
P 3 € l l = R3 12 l/SR 
RE TURN 
ENO 
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SUBROUTINE 

RING A IX) 



COMMON OR( 9) ,K(9) ,ORV (5 ) ,RV ( 5) ,<U (o J , Rkl la ) , K2 (10 I , HR 2 ( 10 ) ,ft 3( 2 >, 
1RR 3 1 2 ) * RA 1 5 ) ,RR4 (5» ,R5 I 1 8 > ,RR5 ( 18 > , R6 Cl 1 I , RR6 I 1 l ) , P l , XK , A A, XII 51, 
2XRI 2 I, E1,E2,N,NB,NBC,M,MB,NSA 
XN=X*N 
>N8=X*NB 
SR = R S HL IX ) 

Z=F 1 *02 

RC=SR»Z 

SRT = RS HLT I X 1 

RCT=-SRT/IRC*RC» 

IFINSA) 1,1,2 

1 CN = CCSIXN) 

SN= SI N I XN) 

CNB = C 0 S I XN B I 
SNB= SINIXNBI 
GO TC 3 

2 CN=> SINI XN1 
SN= COSI XN) 

CNB»= SINIXNB) 

SNB = COSIXNB) 

RC T = -RC T 

3 SS=SN*SN8 
SC = SN*CNB 
CS=CN*SNB 

RAID =PCT*RCT*SS/RC 
PAl 5)<=RCT*CS/PC 
RAI3)*RA(5I/PC 
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♦ DESCRIPTION 
El 


♦ SUBROUTINE RING5IX) ♦ 

4 ■ ' ' 4 

♦ PURPOSE ♦ 

4 CALCULATE AT A GIVEN VALUE OF X THE VALUES OF THE INTEGRANDS ♦ 

♦ OF THE FIFTH SET OF CIRCUMFERENTIAL INTEGRALS OF THE RING 4 

4 4 

4 USAGE ♦ 

4 IT IS USED AS AN ARGUMENT OF THE NUMERICAL INTEGRATION 4 

♦ SUBROUTINE ‘GAUSS* ♦ 

« 4 

IF THE PARAMETERS 4 

7-DISTANCE OF THE SHEAR CENTER OF THE RING FROM THE 4 

MIOOLE SURFACE OF THF SHELL ♦ 

Z-OISTANCE OF THE CENTROtD OF THE RING FROM ITS 4 

SHEAR CENTER ♦ 

THF INPUT VALUE OF X A T WHICH THE INTEGRANDS HAVE 4 

TO BE EVALUATED ♦ 

INTEGRANDS ♦ 

N 4 

NB ♦ 

0 WHEN COMPUTING THE SYMMETRIC MODE SHAPES ♦ 

1 WHEN COMPUTING THE ANTISYMMETRIC MODE SHAPES 4 

♦ 

V 4 

EXCEPT FOR X ALL THE INPUT AND OUTPUT PARAMETERS ARE 4 

COMMUNICATED THROUGH THE COMMON STATEMENT 4 


♦ 

'4 

4 
4 
4 
4 
4 
4 
4 
4 
4 
4 

4 REMARKS 

4 


F2 


R5II! 

N 

NB 

NSA 


♦ FUNCTION SUBROUTINES REQUIRED 
4 R SHt 

4 R ROT 

4 « SHL T 

4 


SUBROUTINE RING5IXI 

COMM CN UK (9) ,R (9) ,UKV (5 ) , «V 15 I ,Kl (8 ),rthl (8 ), *21 10 >, Kh2( 10) ,K 31 2) , 
1RR3I2I .RMS) ,RPM5) . R5I18) , RR5 1 1 8 ) , R6 1 1 1 ),RR6( 11) ,P I , XK , AA. X I ( •>) , 
2XR I 7 ) » E 1,E 2.N.NB.NBC .M.MB.NSA 
XN=X*N 
XNB=X*NB 
SR = R SHL ( X) 

7=E14E2 
RC=SR47 
SRT = R$HLT ( X ) 

RC2 =RC*RC 
RC T =- SR T/RC 2 
BT=RRRT IX ) 

IF (NSA) l ,1 ,2 

1 CN= COSCXN) 

SN= SINIXN) 

CNB = COS(XNB) 

SNB= SINIXNB) 

GC TC ? 

2 CN= SI N I XN) 

SN= COS I XN I 
CNR= SINIXNB) 

SNB = CO SI XNB ) 

RT=-RT 

PCT=-RCT 

’ S S = S N* S N B 
SC= SN*C NB 
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cs=CN*SKe 

U=RCT/IRC2*SR) 

R5(3)=RCT*RCT*SS/(RC*SR » 

R5C2 )=RSC3 »/SR 

R5C 18)=SS*RT*RT/RC 

R5 H )=R5 1 18 I/RC2 

R514)=SS*RT*RCT/RC2 

R5(5)=R5(4)/SR 

R5(6)=SC*RT/RC 

R5l 7)sC S*RT/RC 

R5(ll»=R < 5Cy )/SR 

R5 (10) =R5 (6 )/SR 

R'5( 8)-R5( 61 /RC 2 

R5(9)=R5(7)/RC2 

R5 ( l 2) =R5(8)/SR 

R5( 13)=RS(,9)/SR 

P5 ( 16 >=S C*U 

R5( 17) =CS*U 

R5( 1 4) = R 5( 16) /SR 

R5 ( l 5 ) =R5 (17 )/SR 

RETURN 

ENC 
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♦ USAGE 
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SUBROUTINE RING6IXI 


El 


62 


R 61 I I 
N 

Ne 

NSA 


Z-OISTANCE OF THE SHEAR CENTER OF THE RING FROM THE ♦ 
MIDDLE SURFACE OF THE SHELL 

Z-OISTANCE OF THF CENTROID OF THE RING FROM ITS 
SHEAR CENTER 

THF INPUT VALUE OF X A T WHICH THE l NTEGRANOS HAVE 
TO BE EVALUATED 
INTEGRANDS 
N 
NB 

0 WHEN COMPUTING THE SYMMETRIC MQOE SHAPES 

1 WHEN COMPUTING THE ANTISYMMETRIC MODE SHAPES 


REMARKS 

EXCEPT FOR X ALL THE INPUT AND OUTPUT PARAMETERS ARE 
COMMUNICATED THROUGH THE COMMON STATEMENT 

FUNCTION SUBROUTINES REQUIRED 
R SHL 
RRRT 
RSHLT 


4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 


CCM.MLN OR (9 » ,R19 1 t UKV (5 I, t ,Kkl ( 6i,K2( 10 ) , kh 21 1 0> ,h 31 2 ) , 

1RR3 (2) t R*(5),RR4(5 l,R5( 18 I , RP SI 1 8 > , P. 6( 11 ) ,RR 6(11 ) ,P I , XK , A A , X I ( 51 , 
2XRI2I »E1 »E2 *N»NB»NBC#M,M8«NSA 
XN= X *N 
XN8=X*NB 
SR =R SHL ( XI 
Z= E 1 4£2 
RC=SR4Z 
SRT=RSHLT( XI 
PC2=RC*RC 
RCT.-SRT/RC2 
RT=RRRT( X) 

IF (NSA) If 1« 2 

1 CN* COS ( XN) 

SN= SIN(XN) 

CNB= COS (XNBI 
SNB * SIM XNBI 
GO TO 3 

2 CN= S I N ( XN I 
SN= COS ( XN) 

CNB= SIN(XNB) 

SNB = COS (XNBI 

RT=-RT 

<»CT=-RCT 

3 SS=SN*SNe 
*0 = SN*CNB 
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CS=CN*SN8 
U=RCT/(RC*SR) 
V=RT/RC 2 

P6 ( 1 »*SS*RT *RT/RC2 
R6< 2 > * SS*RT*RCT/RC 
R6I2)=R6( 2»/SR 
R6t4> = SOV 
R6(5»=C S*V 
RM6 >=R6 ( 4 ) /SR 
R6<7»=R6(5)/SR 
R6< ICJ = SC*U 
R6 ( 1 1 )=CS *U 
RM8> =R6<10)/SR 
R6( 9 > *R 61 1D/SR 
RETURN 
END 
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c 
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4 

c 

♦ 

SUBPCUT INE 

GAUSS 

4 

c 

♦ 



4 

c 

♦ 

PURPOSE 


4 

0 

4 

TO EVALUATE THE INTEGRALS OF A SET OF FUNCTIONS OVER THE 

4 


4 

INTERVAL 

LOWL TO UPL 

4 

c 

♦ 



4 

c 

4 

USAGE 


4 

c 

♦ 

CALL GAUSSING, KG, LOWL, UPL, NOFN, PHI ,FX ,SUM ,ANS ,FOFX> 

4 

c 

♦ 

PARAMETER FOFX REQUIRES EXTERNAL STATEMENT 

4 

c 

4 



4 

c 

♦ 

DESCRIPTION 

OF THE PARAMETERS 

4 

c 

4 

NG 

MJMBFR CF POINTS OF THE GAUSSIAN QUADRATURE 

4 

c 

4 

KG 

NUMBER OF SUBINTERVALS 

4 

c 

4 

LOWL 

LOWER LIMIT OF THE INTEGRATION 

4 

r 

4> 

UPL 

UPPER LIMIT OF THE INTEGRATION 

4 


♦ 

NOFN 

NUMBER OF FUNCTIONS TO BF INTEGRATED 

4 

C 

4 

PHI 

AN INTERMEDIET VALUE 

4 

c 

4 

FX 

NAME OF THE VECTOR WHOSE ELEMENTS ARE THE FUNCTIONS 

4 

c 

4 


WHICH ARE TO BE INTEGRATED 

4 

r 

4 

SUM 

WORKING VECTOR 

4 

r 

4 

ANS 

THE NAME OF THE OUTPUT VECTOR WHICH CONTAINS THE 

4 

r 

♦ 


VALUES OF THE INTEGRALS 

4 

r 

4 

FOFX 

THE NAME OF THE EXTERNAL SUBROUTINE USED TO EXPRESS 

4 

C. 

4 


THE FUNCTIONS FX 

4 

c 

4 

TEMP 

WORKING VECTOR WHOSE DIMENSION MUST BF GREATER THAN 

4 

c 

4 


OR EQUAL TO NOFN 

4 

c 

4 

CONST 

VECTOR CONTAINING THE ABSISSAS AND WEIGHT FACTORS 

4 

r 

4 

NLF 

WORKING VECTOR 

4 

c 

4 



4 

c 

4 

REMARKS 


4 

c 

4 

THE 


4 

C 

4 

THIS SUBROUTINE IS EQUIPPEC TO PERFORM 3-, *-, 5-, 6-, 7- , 8- , 9-, 

4 

c 

4 

10-, 16-, 

32-POINT GAUSSIAN QUADRATURE WITH ANY SPECIFIED KG 

4 

c 

4 

NUMBER OF SUBINTERVALS 

4 

c 

4 



4 

c 

4 

subroutines 
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4 

c 

4 

THE EXTERNAL SUBROUTINE FOFX CONTAINING A SET OF FUNCTIONS 

4 

C 

r 

4 
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4 

c 
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SUBROUTINE GAUSS <NG.KG,LOWl,UPL,NCFN,PHl,FX,SUM, ANS, FOFX) 
REAL LOWL 

01 PE NS I ON SUHIl) . ANS ( 1 ) . FX 111 , CONST <1 04 ) , T EMP ( 35 ) ,Nl S C 32 ) 

t 

CATA NLS /-l, -1, 0, *, e, 1*. 20, 28, 36, *6, -1,-1, -1,-1, -1,56, 

* -1,-1, -l, -I, -1,-1, -1, -1,-1, -1,-1, -1,-1, -1,-1, 72 / 

CONST! 1 > =0. 0 

CONST ( 2 ) =0.888888888888889 
COKST13) =0. 77*5966692*2*83 
CONST! *>=0. 555555555555556 
CONST 15 ! =0.33 99 BIO* 3584856 
CONST! 6 1 =0. 6 521*5154 862 5*6 


) 
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CONST (7) =0.86 11 3631 1594 053 
CONST! 8) =0.347854845 137454 
CONST ( 9 1= 0 . 0 

CONST! 10! =0.56 888888 88 88889 
CON STUD =0.538469310105683 
CONST! 121=0.478628670499366 
CONST! 13 » =0 . 906179845938664 
CONST! 141=0.236926885056189 
CONST! 151=0.238619186083197 
CONST! 161=0.467913934572691 
CONST! 171=0.661209386466265 
CONST! 181=0.360761573048139 
CONST! 191=0.932469514203152 
CONST! 201=0.17132449237917 
CONST! 211=0.00000000000000 
CONST! 221=0.417959183673469 
CONST! 231=0.405845151377397 
CONST! 241=0.381030050505119 
CONST! 251=0.741531185595394 
CONST! 261=0.279705391489277 
CONST! 271=0.949107912342759 
CONST! 281=0.12948496616887 
CONST! 291=0.18343464249565 
CONST! 301=0.362683783378362 
CONST! 31 1=0.525532409916329 
CONST 1 321=0 .31 3706645877887 
CONST! 331=0.796666477413627 
CONST! 341=0.222381034453374 
CONST! 351=0.960289856497536 
CONST! 361=0.101228536290376 
CONST! 371=0.000000000000000 
CONST! 381*0.330239355001260 
CONST! 391=0.324253423403809 
CONST! 401=0.312347077040003 
CONST! 411=0.613371432700590 
CONST! 421=0.260610696402935 
CONST! 431=0.836031107326636 
CONST! 441=0.180640160694857 
CONST! 451=0.968160239507626 
CONST! 461=0.081274388361574 
CONST! 471=0.148874338981631 
CONST! 481=0.295524224714753 
CONST! 491=0.433395394129247 
CONST! 501=0.269266719309996 
CONST! 511=0.679409568299024 
CONST! 521 =0.2190863625 15982 
CONST! 531=0.865063366688985 
CONST! 541=0.149451349150581 
CONST! 551=0.973906528517172 
CONST! 561=0.066671344308688 
CONST! 571=0.095012509037637 
CONST! 581=0.189450610455068 
CONST! 591=0.281603550779259 
CONST! 601=0.182603415044924 
CONST! 611=0.458016777657227 
CONST! 621=0.169156519395003 
CONST! 631=0.617876244402644 
CONST! 641=0.149595988816577 
CONST! 651=0.755404408355003 
CONST! 661=0.124628971255534 
CONST! 671=0.865631202387832 
CONST! 681=0.095158511682493 
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CONST! 691=0.944575023073233 
CONST ( 701=0. 062253523938648 
CONST ( 71 )=C. 989400934991650 
CONST ( 721=0.027152459411754 
CONST! 731=0.048307665687738 
CONST! 741=0.096540088514728 
CONST! 751=0.144471961582796 
CONST! 761=0.095638720079275 
CONST! 771=0.239287362252137 
CONST! 781=0.093844399080805 
CONST! 791=0.331868602282128 
' CONST! 801=0.091173878695764 

CONST! 811=0.421351276130635 
CONST! 821=0.087652093004404 
CONST! 831 =0.506899908932229 
CONST! 841=0.083311924226947 
CONST! 851=0.587715757240762 
CONST! 861=0.078193895787070 
CONST! 871=0.663044266930215 
CONST! 881=0.072345794108849 
CONST! 891=0.732182118740290 
CONST! 901=0.065822222776362 
CONST! 911=0.794483795967942 
CONST! 921=0.058684093478536 
CONST! 931=0.849367613732570 
CONST! 941 =0. 050998059262376 
CONST! 951=0.896321155766052 
CONST! 961=0. 042835898022227 
CONST! 971=0.934906075937740 
CONST! 981=0.034273862913021 
CONST! 991=0.964762255587506 
CONST! 1001=0.025392065305262 
CONST 1101 1 = 0.985611511545268 
C0NST11021 =0.016274394730906 
CONST! 1031=0.997263861845482 
CCNST 1104 1=0.007018610009470 
1001 FORMAT I31HIG6USS TAN INTEGRATION OF ORDFR 15, / 

* 17H01S NOT AVAILABLE / 

* 1 7 HO PROGRAM STOPPED. 1 
C 

C FROM ORCER SPECIFIED, FIND START LOCATION IN OATA. 
00 20 N0R0=1 ,32 
IF (NG .EQ. NORO) GO TO 25 
20 CONTINUE 
99 1.RITE (6,10011 NG 
CALL EXIT 
25 LS = NLS (NORO I 

IF (LS .LT. 0) GO TO 99 
FKG= FLOATING! 

00 l N= 1 , NOFN 
1 ANS ! N1 =0.0 
A1=L0WL 
A2 = A1 
B=UPL 

0X= < B-Al l/FKG 
CO 2 I =1 , KG 
A 1 = A2 
e=Ai40x 
A2 =8 

At P HA = C. 5E 0*0 X 
BET A= AL PHA*A l 
00 3 N=l ,NOFN 
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3 SUMI.N >=0.0 

0C 4 J=1,NG,2 
LCC = LS ♦ J 

PYKA = CONST! LOCI * ALPHA 
PH 1 = BET A-PYK A 
CALL FOFX ( PH I I 
CO 5 N= 1 * NOFN 

5 TEP P 1 N I =FX ( N I 

IF(PYKA) 7,6,7 

6 CO 8 N= 1# NOFN 
8 TEMPINI=0.0 

GO TO 9 

7 PH!=PYKA*BFTA 
CALL FOFX (PHI) 

S 00 10 N= l ,NOF N 
L0CP1 = LS ♦ J*1 

10 SUM ( NJ = SUM(N) ♦ (TEMP1NI ♦ FX(NI) 

4 CONTINUE 

CC 11 K= 1 » NCFN 

11 ANSINI =ANS(Nl4SUM(NJ*AL PHA 
2 CONTINUE 

RETURN 

ENO 


CONST I LOCP 1 1 
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O <-> W VJ o u 


SUBROUTINE EIGENC I AK , AK , MN3 , KRRR, V ECR, L C,XXXX, r , MC , Z , EVR , EV I , I NOI 
1C I 

Cl RE NS ION AM (KRRR, 11 , AKIKRRR, 1 1 , VECR < KRRR , 1 1 ,LC ( II ,E VRC 1 1 .XXXXIl) , 

mu ,mc( 1) ,zm , iNDicn i.evi m 


the purpose of this subroutine is to arrange the matricies into 

the FORM REQUIRED BY SUBROUT1 NE EIGENP 


CALL RRAY( 2, MN 3, MN3, KRRR, KRRR, XX XX, AM) 

CALL RRAY (2,MN3,MN3, KRRR, KRRR, Y,AK) 

CALL INV(XXXX,MN3,LC,MC) 

CALL MPRO(XXXX,Y,Z,HN3,MN3,MN3> 

CALL RRAY (l, MN3.MN3, KRRR, KRRR, Z,AK) 

00 2C05 I J=l , MN3 
EVR< IJ 1 = 0.000 
FV I { I J ) =0 .000 
00 2C05 KJ=1 ,MN3 
AM( IJ.KJ 1=0. COO 
20CS VFCR(IJ,KJ1=0.000 
U=4B. 000 

CALL E IGENP(MN3,KRRR .AK.U.E VR.E VI , VECR, AM, I NOI C) 

RETURN 

PND 
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SUBROUTINE E IGENPI N,NM, A ,T ,E VR ,EV I ,VECR , VEC I , INOIC I 
C THIS SUBROUTINE TAKEN FROM COMMUNICATIONS OF ACM VOL. II NO. 12, 

C DEC. 68 

INTEGER I ,1 VEC,J,K,Kl,KON,L,Ll,M,N,NM 
DIMENSION A (NM, 1 ) , VECR I NM, 1 >, VEC I (NM, 1 > , 

IFVRINM) ,EVI (NM) , I NO I C ( NM ) 

DIMENSION I WORK! 120) , LOCAL I 120) , PRFACTM20) 

1 , SUBDI AI 120 ) ,W0RK1 1 120 ) , WORK 21 120 ) , WORK! 120) 

C 

C THIS SUBROUTINE FINOS ALL THE EIGENVALUES ANC THE 
C EIGENVECTORS OF A REAL GENERAL MATRIX OF ORDER N. 

C 

C FIRST IN THE SUBROUTINE SCALE THE MATRIX IS SCALED SC THAT 
C THE CORRESPONDING ROMS ANO COLUMNS ARE APPROXIMATELY 
C BALANCED ANO THEN THE MATRIX IS NORMALISED SC THAT THE 
C VALUE OF THE EUCLIOIAN NORM OF THE MATRIX IS EQUAL TC CNE. 

C 

C THE EIGENVALUES ARE COMPUTEO 8Y THE 03 DCUBL E-ST EP METHOD 
C .IN THE SUBROUTINE HESOR . 

C THE EIGENVECTORS ARE COMPUTED BY INVERSE ITERATION IN 
C THE SUBROUTINE REALVE, FOR the REAL EIGENVALUES, CR IN THE 
C SUBROUTINE COMPVE,FOR THE COMPLEX EIGENVALUES. 

C 

C THE ELEMENTS OF THE MATRIX ARE TO BE STORED IN THE FIRST N 
C RCWS- ANC COLUMNS OF THE TWO DIMENSIONAL ARRAY A. THF 
C CRIGINAL MATRIX IS OESTRCYEO BY THE SUBROUTINE. 

C N IS THE OROER OF THE MATRIX. 

CAM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL 
C ARRAYS A , VECR , VEC I AND THE DIMENSION OF THE CNF 
C DIMENSIONAL ARRAYS EVR.EVI AND INOIC. THE RE F CRB THE 
C CALLING PROGRAM SHOULC CONTAIN THE FOLLOWING. DECLARATION 
C DIMENSION AINM.NNI ,VECR(NM,NN) ,VECI INM.NN), 

C 1 EVR ( NM ),EVI(NM), INDIC(NM) 

C WHERE NM ANO NN ARE ANY NUMBERS EQUAL TO CR GREATER THAN N 
C THE UPPER LIMIT FOR NM IS EQUAL TO 100 BUT MAY BE 
C INCREASED TO THE VALUE MAX BY REPLACING THE DIMENSION 
C STATEMENT 

C DIMENSION IWORKI 100), LPCAL( 100) , . .. ,W0PKI 100) 

C IN THE SUBROUTINE EIGENP WITH 

C DIMENSION IWORK(MAX) .LOCAL(MAX) .... ,WORK(MAX) 

C NM AND NN ARE OF COURSE BOUNDED BY THE SIZE OE THE STORE. 

C 

C THE REAL PARAMETER T MUST BE SET EQUAL TO THE NUMBER CF / 

C BINARY DIGITS IN THE MANTISSA PF A DOUBLE PRECISION 
C FLOATING-POINT NUMBFR. 

C 

C THE REAL PARTS OF THE N COMPUTEO EIGENVALUES WILL BE FOUND 
C IN THE FIRST N PLACES OF THE ARRAY EVR AND THE IMAGINARY 
C PARTS IN THE FIRST N PLACES OF THE ARRAY EVI . 

C. THE REAL COMPONENTS OF THF NORMALISED EIGENVECTOR I 
C <1=1 ,2, . . . ,N( CORRESPONDING TO THE EIGENVALUE STORED IN 
C EVR < I ) AND EVI(I) WILL BE FOUND IN THE FIRST N PLACES OE 
C THE COLUMN I OF THE TWC DIMENSION AL ARRAY VECR AND THE 
C IMAGINARY COMPONENTS IN THE FIRST N PLACES OF THE COLUMN I 
C TF THE TWO DIMENSIONAL ARRAY VEC I. 
f 

C THE REAL EIGENVECTOR IS NORMALISED SO THAT THE SUM OF THE 
C SCUARES OF THE COMPONENTS IS EQUAL TO ONE. 

C THE COMPLEX EIGENVECTOR IS NORMALISED SO THAT THE 
C COMPONENT WITH THE LARGEST VALUE IN MODULUS HAS ITS REAL 
C PART EQUAL TO ONE ANC THE IMAGINARY PART EQUAL TO ZERO. 

C 
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n n o n o rv. 


C THE ARR AY ! NCIC INDICATES THE SUCCESS OF THE SUBROUTINE 

EIGENVECTOR I 
NOT FOUND 
NCT FOUND 
FOUND 


IFIN.NE. 1IGO TO ' 1 
EVP III = Ml, II 
EVI ( II = 0. 00 
VECR ( 1, 11 = l.OC 

VECM1.1 I = 0.00 s 

INDICIl) = 2 
CO TO 25 
C 

1 CALL SCALE! N.NM.A.VECI .PRFACT.ENORM) 

C THE COMPUTATION OF THE EIGENVALUES OF THE NORMALISED 
C MATRIX. \ 

EX = EXPI-T*AL0G!2. COI » 

CALL HESORIN.NM, A.VECI, EVR, EVI, SUBDI A , IND IC ♦ EP S.E X I 

C 

C THE POSSIBLE DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX 
C INTO THE SUeMATRICES OF LOWER ORDER IS INDICATED IN THE 
C. ARRAY LOCAL. THE DECCMP0S1T I ON OCCURS WHEN SCME 
C SUBDIAGONAL ELEMENTS ARE IN MODULUS LESS THAN A SMALL 
C POSITIVE NUMBER FPS CEF INED IN THE SUBROUTINE HESQR . THE 
C AMOUNT OF WORK IN THE EIGENVECTOR PROBLEM MAY BE 
C C1MINISFED IN THIS WAY. 

J = N 
I = l 

LOCAL (1 I = 1 
IF! J.EO.l IGO TO 4 

2 IF! ABS! SU8DIA! J-U'.GT.EPSIGC TO 3 
I = !♦! 

LOC AL ( I ) =0 

3 J = J-l 

LOCAL! I ) = LOCAL! Iltl 
IFIJ.NE.IIGO TO 2 

C 

c THE EIGENVECTOR PROBLEM. 

4 K = 1 
KON = C 

L = LOCAL! I 1 
M * N 

DO 10 1 = 1, N 
IVEC = N- I *1 
IFd.LE.LIGO TO 5 
K = K ♦ 1 
M = N-L 

L = L«- LOCAL ! KI 

5 IF! INOIC! IVEC I .EQ. CIGO TO 10 
IFtEVH IV FCI.NE, 0.001 GO TO 8 

C 

C TRANSFER OF AN UPPER-HESSENBERG MATRIX OF THE ORDER M FROM 
C THE ARRAYS VECI ANO SUBCIA INTO THE ARRAY A. 

DC 7 Kl=l,M 
DO 6 Ll=Kl,H 

6 A I Kl , LI I = VECIIKltLlI 
IF! Ki.EQ. 1 IGO TO 7 

A ! K l » K 1- l I = SUBDI A ! K 1- 1 ) 

7 CONTINUE 


EIGENP AS 

FOLLOWS 


VALUE 

OF INOIC! 1 1 

EIGENVALUE t 


0 

NCT FOUND 


1 

FOUND 


2 

FOUND 
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C THE COMPUTATION OF THE REAL EIGENVECTOR I VEC OF THE UPPER- 
C HESSEN8ERG MATRIX CORRESPONDING TO THE REAL EIGENVALUE 
C EVRIIVECI. 

CALL REAL VEIN* NM,M, IV EC, A, VECR,E VR ,E VI .IWORK, 

1 WORK, INDIC.EPS, EX) 

GO TO 10 
C 

C THE COMPUTATION CF THE COMPLEX EIGENVECTOR IVEC OF THE 
C UPPER-HESSENBERG MATRIX COR RE SPONOt NG TO THE COMPLEX 
C EIGENVALUE EVRIIVECI ♦ 1*EVI(1VECI. IF THE VALUE OF KON IS 
C NOT EQUAL TO ZERO THEN THIS COMPLEX EIGENVECTOR HAS 
C ALREADY BEEN FOUNO FROM ITS CONJUGATE. 
f IF I KON. NE . 0 IGO TO 9 
KON = 1 

CALL CCMPV E1N«NM,m, IV EC, A, V ECR, V EC I , EVR , EV I , INO IC , 

1 I WORK « SU80 1 A , WOR Kl ,M0RK2 , WORK, EPS, EX) 

GO TO 10 
9 KCN = 0 
1C CONTINUE 
C 

C THE RECONSTRUCT ION OF THE MATRIX USEO IN THE REDUCTION OF 
C MATRIX A TO AN UPf>E»-HE SSENBERG FORM BY HCUSEHOLCFR METHOC 
CO 12 1=1, N 
OC 11 J=I ,N 
A I I , J I = 0. 00 

11 AIJ, I ) = 0.00 

12 All ,1 ) = 1.00 
IFIN.LE.2IG0 TO 15 
M = N-2 

00 14 K = 1 , M 
L = KM 
DC 14 J=2,N 
01 = 0.00 
00 13 t=L,N 

02 = VECII I , K ) 

13 01 = DU- 029AI J,I I 

00 14 I=L,N 

14 A I J , 1 1 = AIJ, l I-VECII I,KI*Ol 
C 

C THE COMPUTATION OF THE EIGENVECTORS OF THE ORIGINAL NON- 
C SCALED MATRIX. 

15 KON = 1 

CO 24 1=1, N 
L * 0 

IFIEVII I I .EQ.O.OOIGO TO 16 
L = l 

I F (KON . EQ.O I GO TO 16 
KON = 0 

GO TO 24 

16 00 18 J=1,N 

01 = 0.00 
02 = 0.00 
OC l? K=l,N 

03 = A I J , K I 

01 = 01 *03 *V ECR IK, II 
I F I L. EQ.O) GO TO 17 

02 = D2*03*VECR(K,l-l) 

17 CONTINUE 

WORK I J I = 0 1/PRF ACT I J I 
IFIL.EQ.OIGO TO 18 
SUBOIAI JI=02/PRFACT( Jl 
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18 


CONTINUE 


C TEE NORMALIZATION OF THE EIGENVECTORS ANO THE COMPUTATION 
C OF THE EIGENVALUES OF THE ORIGINAL NON-NORM AH SEO MATRIX. 
IFIL.E0.1IG0 TO 21 
01 = 0.00 
OC 19 M=1 ,N 
,H1 = WORklMH 

19 ' Ol = 0UW1*W1 ‘ 

01 = SORT (01 I 

00 20 M=1 ,N 

VEC KM, I) = 0.00 

20 VECRIH.n = WCRMMI/D1 
E VR ( I ) = E VR ( I 1 *E NORM 

GC TO 24 
C 

21 KON = 1 

EVR I I ) = EVR(!)*ENOBM 
EVRII-1I = FVR(I) 

C V I ( I I = EVK I) *ENORM 
FVI (.1-1 I =-EV I ( I I 
R = C.00 
CO 22 J= 1 » N 
W1 = WORK ( J ) 

W2 = SUBOIA(J) 

R l = W1*W1 ♦ W2*N2 
IF(R.GE.R1 IGO TO 22 
R = R 1 
L = J 

22 CONTINUE 
03 = WORK! L I 

n 1 = SUBO I A (L I 
00 23 J«1,N 

01 = WORK ( J I 

02 = SUBOIA(J) 

VECRIJ.I) = (01*03*-D2*R1I/R 
VECKJ.tl = ( 02*03-0 l *R 1 I /R 
VECRI J.I-l I = VECRIJ, I ) 

23 VEC I ( J, I - 1 1 *- VEC I ( J , 1 1 

24 CONT INUE 
C 

25 RETURN 
ENC 
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SUBROUTINE SCALE < N ,NM ,A ,H , PRF ACT , ENCPM) 

INTEGER I» J * ITER ,N, NCOUNT » NM 

Cl PENS ION A(NM,1 I ,HINM, H.PRFACTINM) 

C 

C IMS SUBROUTINE STORES THE MATRIX OF THE ORDER N FROM THE 
C ARRAY A INTO THE ARRAY H. AFTERWARD THE MATRIX IN THE 
C ARRAY A IS SCALED SO THAT THE QUOTIENT OF THE ABSOLUTE SUM 
C CF THE OFF-Cl AGONAL ELEMENTS OF COLUMN I AND THE ABSOLUTE 
C SUM CF THE OFF-DIAGONAL ELEMENTS OF ROW I LIES WITHIN THE 
f. VALUES OF BOUNOi ANO BOUND2. 

C THE COMFCNENT I OF THE EIGENVECTOR OBTAINEC BY USING THE 
C. SCALED MATRIX MUST BE DIVIDED BY THE VALUE FCUNO IN THE 
C PRFACTCII OF THE ARRAY PR FAC T . IN THIS WAY THE EIGENVECTOR 
C OF THE NON-SCALEO MATRIX IS CBTAINED. 

C 

C AFTER THE MATRIX IS SCALED IT IS NORMALISED SO THAT THE 
C VALUE OF THE EUCLIDIAN NCRM IS EQUAL TO ONE. 

C IF THE PROCESS OF SCALING WAS NOT SUCCESSFUL THE ORIGINAL 
C MATRIX FROM THF ARRAY H WOULD BE STORED BACK INTO A ANO 
C THE EIGFNPRCBLEM WOULD BE SOLVED BY USING THIS MATRIX. 

0 NM CEFINES THE FIRST DIMENSION OF THE ARRAYS A AND H. NM 
C MUST BE GREATER OR EOUAL TO N. 

C THE EIGENVALUES OF THE NORMALISED MATRIX MUST BE 
C MULTIPLIED BY THE SCALAR ENORM IN ORDER THAT THEY BECOME 
C THE EIGENVALUES OF THE NCN-NCRMAL IS ED MATRIX. 

C 

CO 2 1=1, N 
DC 1 J=l , N 
l HII.JI =A f I , J) 

Z PR FACT 111= 1.00 
BOUNDl = 0.7500 
B0UNC2 = 1.3300 
ITER = 0 
3 NCOUNT = 0 
CO 9 1=1, N 

COLUMN = 0.00 
ROW = C. OC 
CO A J=l,N 

IF I I .EQ. J I GO TO A 

COLUMN = COLUMN ABSIAIJ, I)> 

ROW = ROW ♦ ABSI A( I.JII 
A CONTINUE 

I F ( COLUMN.EQ. 0. 00) GC TC 5 
IF I ROW . EQ .0 • 00 I GO TO 5 
Q = COLUMN/RCW 
IFIQ.LT. BOUNDl >G0 TO 6 
IF (Q.GT , 80UNC2 JGO TO 6 

5 NCOUNT = NCOUNT ♦ l 
GO TO 8 

6 FACTOR = SORT (Q ) 

DO 7 J=l,N 

IF( I .EQ .J IGO TO 7 
All ,J> = A ( I , J ) *F ACT OR 
A( J , I I = A< J.II/FACTCR 

7 CONTINUE 

PRF ACT ( 1 1 = PRFACT(M*F ACTOR 

8 CONTINUE 
ITER = IT ER*1 

IP I I TER.GT. 301 GO TO 11 
IF! NCOUNT. LT.NIGO TO 3 
C 

FNORM = 0.00 
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CO 9 1*1, N 
DC 9 J=l .N 
0=d< f , Jl 

9 FNORM = FNORM«-Q*Q 

FNORM = SORT I FNORM) 

00 1C 1=1 .N 
OC 10 J=1'*’N 

lC J i > Ad . J» = A C l . JI7FN0RM 
FNORM ' J s -FNORM . 

no re n 

11 00 12 1 = 1, N 

PRFACT ( I )=1 .00 
00 12 J=1 , N 

12 AII.J) = H(1 ,J) 
FNCRM = 1 .00 

12 RETURN 
END 
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SUB R CUT I NE HES QR ( N, NM, A , H, EVR , EV I , SUBO I A • INO IC.EP S, EXI 
INTEGER I, J,K,L,M,MAXST,Ml ,N,KM,NS 

Cf REAS ION A INN, 1 1, H(NM, 1 1, E VR (NM I ,E VI I NM I .SUBOIAINM) 
DIMENSION INCICINMJ 

C 

C T> IS SUBROUTINE FINDS ALL THE EIGENVALUES OF A REAL 
C GENERAL MATRIX. THE ORIGINAL MATRIX A OF ORDER N IS 
C REDUCED TO THE UPPFR-HE SSEN8ERG FORM H BY MEANS CF 
C SIMILARITY TRANS FORMAT IONS! HOUSEHOLDER METHOD). THE MATRIX 
C H IS PRESERVED IN THE UPPER HALF OF THE ARRAY H AND IN THE 
C ARRAY SUBOIA. THE SPECIAL VECTORS USED IN THE DEFINITION. 

C CF THE HCUSEHOLOER TRANSFORMATION MATRICES APE STORED IN 
C THE LOWER PART OF THE ARRAY H. 

C NM IS THE FIRST DIMENSION OF THE ARRAYS A AND H. NN MUST 
0 BE EQUAL TO OR GREATER THAN N. 

C THE REAL PARTS OF THE N EIGENVALUES WILL BE FCUNO IN THE 
C FIRST N PLACES OE THE ARRAY EVR, AND 

C THE IMAGINARY PARTS IN THE FIRST N PLACES. OF THE ARRAY FV I 
e the ARRAY INOIC INDICATES THE SUCCESS OF THE ROUTINE AS 
C F C'l'L'CWS 

C VALUE OF INOICIII EIGENVALUE I 

C 0 NOT FOUND 

C 1 FOUND 

C EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS 
C ZERO IN THE PROGRAM. EPS - (EUCLIDIAN NORM OF H»*EX .WHERE 
C EX = 2**(-T>. T IS THE NUMBER OF BINARY OIGITS IN THE 
C MANTISSA OF A FLOATING POINT NUMBER. 

C 

C 

r 

C RECUCTION OF THE MATRIX A TO AN UPPER-HESSENBERG FORM H.. 

C THERE ARE N-2 STEPS. 

IFIN-21 1A, 1,2 

1 SUBOIAU ) = A ( 2, I) 

GO T.C 14 

2 M = N-2 

re 12 K=l,M 
L = K*1 
S = 0.00 
DC 3 I=L,N 

HU , K) - A ( I , K ) 

3 S = S* AB S( A ( I , K ) I 
IFCS.NE. A8S(A|K»l*K> I ICO TO « 

SUBDIA(K) = A ( *♦! ,K) 

F ( K ♦ 1 , K I = 0.00 

GC TO 12 

A SR 2 = 0.00 
. DO 5 I=L,N 

SR = A( I , K I 
SR * SR/S 
A ( I,K I = SR 

5 SR2 = SR2*SR*SR 
SR =■ SQRTISR2I 
IF. (A(L.K).LT .0.00) GO TO 6 
SR = -SR 

t SR 2 * SR 2- SR* AIL, K I 
ACL.KI = AIL.KI-SR 
H(L,K) = H ( L • K I -S R*S 
SUBDIAUO = SR*S 
X = S* SORT (SR2) 

DC T I«L,N 

HI I ,K I = HU ,KI ft. 
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7 SUBC!A( II = AI !,K I/SR2 

C PREMLLTIPLICATION BY THE MATRIX PR. 

CO 9 J=L r N 
SR‘ = 0.00 
00 8 I =L »N 

8 SR = SR«-Al t,K »*AI I, J I 

00 9 1=1. N 

9 A ( l ,'J ) = A I I , J I- SUBD I A ( 1 1 • SR 
C PC ST PUL T I PL I CAT I ON BY THE MATRIX PR. 

DO 11 J=1 ,N 
SR=0 .00 
00 10 I = L » N 


1C 

SR = SR*AI J , I I *A1 I ,KI 



OO'li I = L,N 


11 

AIJ.ll = A I J , I 1-SUBOI A1 1 |*SR 


12 

CONTINUE 



CO 13 K=1,M 


13 

A I K*1 ,K1 = SUBOIAIK l 



C TRANSFER OF THE UPPER HALF OF THE MATRIX A INTO THE 
C ARRAY H ANC THE CALCULATION OF THE SMALL POSITIVE NUMBER 
C EPS. 

SUBC1AIN-1) = AIN.N-ll 

14 EPS = 0.00 
no 15 K=1 ,N 

INC IC I K I = 0 

l F (K . NE.N 1 EPS = EPS ♦ SUBD IAI K I *SUBO I A IK I 
00 15 I=K,N 

MK, I ) = AIK, I ) 

WZ = AIK, I) 

15 EPS = EPS ♦ U2*U2 

EPS = EX* SORT! EPS) 

THE CR ITERATIVE PROCESS. THE UPPER-HES S ENB E PG MATRIX H IS 
RECUCED TO THE UPPER-MOO I F I EO TRIANGULAR FORM. 

DETERMINATION OF THE SHIFT CF ORIGIN FOR THE FIRST STEP OF 
THE OR ITERATIVE PROCESS. 

SHIFT = AIN.N-ll 
IFIN.LE.2I SHIFT = 0.00 
IF! AIN.N ) .NE .O.OOISHIFT = 0.00 
IF! AIN-l ,N» .NE.0.001SHIFT = 0.00 
!F< A1N-1.N-1I.NE. O.OOISHIFT = 0.00 
P = N 
NS = 0 

MAXST = N* 1 0 

TESTING IF THE UPPER HALF OF THE MATRIX IS EQUAL TO ZERO. 
IF IT IS EOUAL TO ZERO THE QR PROCESS IS ACT NECESSARY. 

CC 16 1=2, N 
DO 16 K = I ,N 

I F I A I I-l.Kl.NE.O. 001 GO TO 18 

16 CONTINUE 

00 17 I =1 ,N 
INOICI 11= l 
EVRIII = All, II 

17 E VII 1 1 = 0.00 
GO TO 37 

START THE PAIN LCCP CF THE QR PROCESS. 

18 K=M- 1 
P1=K 

1 = K 
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C FIND ANY DECOMPOSITIONS OF THE MATRIX. 

C JUMP TO 34 IF THE LAST SL8MATRI X OF THE OFCOPCS IT ION IS 
C 01 TM= ORDER ONE. 

C JURE TC 35 IF THE LAST SUBMATRIX OF THE DECOMPOSITION IS 
C OF THE ORDER TWO. 

IFIK 137,34, 19 

19 IF! A8SIAI M, K) I .LE.EPS I GO TO 34 
IFIM-2.EO.OIGO TO 35 

20 I = 1-1 

IF( ABSI AIK, I) ) .LE.EPS IGO TO 2t 
K = I 

IFIK .GT .1 IGO TO 20 

21 IFIK. EO. Ml IGO TO 35 

C TRANSFORMATION OF THE MATRIX OF THE ORDER GREATER THAN TWC 
S = AIM.MUAIMl.Ml I+SHIFT 

SR = AIM,MI*AIM1 ,M1 l-AIM.Ml |*AIM1,M| ♦ 0 . 2500 *SH I FT »SH FT 
AIK ♦2*K I = 0.00 

C CALCULATE X1,Y1,Z1,F0R THE SUBMATRIX OBTAINEC BY THE 
C DECOMPOSITION. 

X = AIK,K)*|A|K,Kl-SUAIK,K*ll*AIKM»KI«-SR 
Y = AIK»1 ,K»*IA{K,KI+AIK*l,K*ll-S I 
R = ABSI XI* ABSI Yl 
IFIR.EQ.O.OOISHIFT = AIM.M-ll 
IFIR.EO.O.OOIGO TO 21 
2 = AI K*2,K* 1I*AIK* l ,KI 
SHIFT = 0.00 
NS = NS*l 

THE LOOP FOR ONE STEP OF THE OR PROCESS. 

DO 33 I = K , Ml 

IF! I.EO.K IGO TO 22 
C CALCULATE XR.YR.ZR. 

X = AI 1,1-11 

Y a AI 1*1, I-1I 

l - 0.00 

I F I U2.GT.MIG0 TO 22 
Z = AI 1*2, 1-1 I 

22 S«2 = ABSIX1* ABS (Y I ♦ ABS1Z1 
t F I SR 2.EQ. 0. OOIGO TO 23 

X = X/SR2 

Y = Y/SR2 
Z = Z / SR 2 

23 S = SORT |X*X ♦ Y*Y *Z*Z I 
IF! X.LT.O. OOIGO TO 24 

S = -S 

24 IFII .EQ.KIGC TC 25 
At I ,1- l) = S*SR2 

25 IF ISR2.NE.0. OOIGO TO 26 
l F 1 1 ♦ 3. GT. Ml GO TC 33 

GO TO 32 

26 SR = 1.00-X/S 
S = X-S 

X = Y/S 

Y = Z/S 

C PREML'LTI PLICATION BY THE MATRIX PR. 

00 28 J*l,M 

S = A 1 1 , J I ♦ A I I *1 , J )*X 
IFI U2.GT.MIG0 TO 2 7 
S = S*At U2, J l*Y 

27 S * S*SR 

AI I , JI = All ,J)-S 
AIU1.JI » AII*l,JI- S*X 
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IF! I+2.GT.MIGC TO 28 
AI I*2,JI = A( I*2,J)-S*Y 
28 CONTINUE 

C POSTMULTIPLICATICN BY THE MATRIX PR. 

L = I ♦ 2 • 

IFU.LT.M1 I GO TO 29 

v . L - M ■ v., , 

2.9 .00 - 31 J*K *L : .., t. . • 

S = AIJ, II+AIJ, t+l)*X 
IF I l*2.GT.M)G0 TO 3C 
S = S ♦ A ( J • I*21*Y 

30 S = S*SR 
A( J,II = A ( J , I I - S 
AIJ»I*1I=AIJ, I ♦ 1 I - S * X 
IF 1 I *2. GT . H) GO TO 31 
AIJ, I*2I=AI J,T*2I-S*Y 

31 CONTINUE 
IFIH-3.GT.MIG0 TO 33 
S = -A( 1*3, I *-2 I ♦ Y* SP 

32 • A ( 1*3 * I ) = S 
AI 1*3 « l ♦ l I = S*X 
A( I+3,I*2) = S*Y ♦ AI 1*3,1*21 

33 CONTINUE 

IF( NS .GT .MAXST IGO TO 3T 
GO TO 18 

COMPUTE THE LAST EIGENVALUE. 

34 EVRIMI = AIM, Ml 
EVIIMI = 0.00 
INCIC(M) = l 
M = K 
GO TO 10 

C COMPUTE THE EIGENVALUES CF THE LAST 2X2 MATRIX OBTAINED BY 
C THE DECOMPOSITION. 

35 P = 0.500*1 AIK.K l+AIM.M I I 
S = 0.500*1 AIM, M)-A(K,K> I 
S = S*S ♦ AI K ,M) *A< M,K> 

INC 1C IK I = 1 

I NO IC I M> = 1 

IF! S.LT.O.OOIGO TO 36 

T = SORT IS I 

FVRIKI = R-T 

FVRIMI = R ♦ T 

FVI IK) = o.no 

F VI I Ml = 0.00 

M = M-2 

GC TC 18 

36 T = SORTI-SI 
EVRIK1 = R 
FVIIK) = T 
FVRIMI = R 

EV 1 1 M I = -T 
H = M-2 

GC TC 18 
C 

37 RETURN 
s FNO 
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SU8R0UTINF R€ALVE(N,NM, M, l VEC . A , VEC R . E VR ,E V! , 

1 1 WORK 'WORK tl NO! C, EPS, EX I 
INTEGER I, I VEC,! TER , J,K ,L ,M,N,NM,NS 
0! PENS ION A (NM , 1 ),VECR(NH, 11, EVRINM ) 

DIMENSION EVIINM) , I WORKING) ,WORK ( NM) , INDICINM) 

C 

C THIS SUBROUTINE FINOS THE REAL EIGENVECTOR OF THE REAL 
C UPPER-HESSENBERG MATRIX IN THE ARRAY A, CORRESPONDING TO 
C THE REAL EIGENVALUE STORED IN EVRIIVEC). THE INVERSE 
C ITERATION METHOD IS USEO. 

C NOTE THE MATRIX IN A I S DESTROYED BY THE SUBROUTINE. 

C N IS THE CROER OF THE UPPEP.-HESSENBERG MATRIX. 

C NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL 
C ARRAYS A ANO VECR . NM MUST BE EQUAL TO OR GREATER THAN N. 
C, M IS THE CROER OF THE SUBMATRIX OBTAINED BY A SUITABLE 
C DECOMPOSITION OF THE UPPER-HESSENBERG MATpiX IF SOME 
C SU8C I AGONAL ELEMENTS ARE EQUAL TO ZERO. THE VALUE OF M IS 
C CHOSEN SO THAT THE LAST N-M COMPONENTS OF THE EIGENVECTOR 
C ARE ZERO. 

C IVFC GIVES THE POSITION OF THE EIGENVALUE IN THE ARRAY E VR 
C FOR WHICH THE CORRESPONDING EIGENVECTOR IS COMPUTED. 

C THE ARRAY EVI WOULD CONTAIN THE IMAGINARY PARTS CF THE N 
C EIGENVALUES IF THEY EXISTED. 

C 

;C THE M COMPONENTS OF THE COMPUTED REAL EIGENVECTOR WILL BE 
£ FCUNC IN THE FIRST M PLACES OF THE COLUMN IVFC OF THE TWO 
C DIMENSIONAL ARRAY VECR. 

C 

C IWCRK ANO WORK ARE THE WORKING STORES USEO DURING T HE 
C GAUSSIAN ELIMINATION ANO BACKSUBST ITUT ION PROCESS. 

C THE ARRAY INDIC INDICATES THE SUCCESS DF THE ROUTINE AS 
C FCILCWS 

C VALUE OF INDICm EIGENVECTOR I 

C l NOT FOUND 

C 2 FOUND 

C EPS ISA SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS 
C ZERO IN THE PROGRAM. EPS = (EUCLIDIAN NORM OF 4)*EX, WHERE 
C EX = 2** ( -T J . T IS THE NUMBER OF BINARY CIGITS IN THE 
C MANTISSA OF A FLOATING POINT NUMBER. 

VEC R < 1 » IV EC ) =» 1.00 
IF(M.EQ.lJGa TO 24 

C SMALL PERTURBATION OF EQUAL EIGENVALUES TO OBTAIN' A FULL 
C SET OF EIGENVECTORS. 

EVALUE = EVRUVECI 
I F ( IVEC .EQ.MIGO TO 2 
K = IVEC*1 
R = C.00 
CO l !=K,M 

I F ( F VALUE . NE. EVR ( 1 1 ICO TO l 
I F ( E V I ( II .NE.O.DOIGO TQ l 
R * RO.OO 

1 CONTINUE 

EVALUE = EVALUE+R *E X 

2 CC 3 K=i,M 

3 A (X , K I = A( K.KI-EVALUE 
C 

C GAUSSIAN ELIMINATION OF THE UPPER-HESSENBERG MATRIX A. ALL 
C ROW INTERCHANGES ARE INDICATED IN THE ARRAY IWORK. ALL THE 
C MULTIPLIERS ARF STORED AS T HF SUBOIAGONAL ELEMENTS OF A. 

K = M-l 
CO 3 1=1, K 
L * I ♦ l 
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I WORK ( I) = 0 

IFtAtI*l, II.NE.O.OOIGO TO A 
IF (At I ,I» .NE.O.OOIGO TO 8 
At 1, I » = EPS 
GC TO 8 

4 I F ( ABStAt I.IM.GE. ABStAt \+l ,11) )GO TC 6 
IWCRKtn = l 

DC 5 J = I ,M 
R * A(!,J) 

At t,JI = At 1*1, j I 

5 A ( I M , J > = R 

6 R = -At IM, t»/A( l , I) 

At 1*1, I) = R 

DC 7 J=L,M 

7 At I*1,J I = At t»l,J)+R*At I ,JI 

8 CONTINUE 

IFtAtM, Ml. NE.O.OOIGO TO 9 
AIM,M) = EPS 

C 

C THE VECTOR tl,l,...,l) IS STORED IN THE PLACE CF THE RIGHT 
C HANC SIOE COLUMN VECTOR. 

9 00 11 I =1 ,N 

IFtl.GT.M) GO TO 10 • . 

WCRKtII = 1.00 
GC TO 11 

1C WORK t II = C.00 

11 CONTINUE 
C 

C THE INVERSE ITERATION IS PERFORMED ON THE MATRIX UNTIL THE 
f. INFINITE NORM OF THE RIGHT-HAND SIDE VFCTOR IS GREATER 
C THAN THE BOUND OEFINEO AS 0.01/tN*EXI. 

BOUNO = C.OlOO/f EX * FL0ATIN1I 
NS * 0 
ITER = 1 

THE PACKSUBSTITUTICN. 

12 R = C.CO 
CO 15 1=1, M 

J = M- I *1 
S = WORK ( J ) 

IFtJ.EO.MIGO TO 14 
L = J*1 
00 13 K=L,m 
SR = WORK t K I 

13 S = S - SR * A t J , K I 

1 4 WORK! J) = S/A t J , J ) 

T = ABS t WORK t J I ) 

IFIR.GE.T IGC TC 15 
R = T 

15 CONTINUE 

C THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR FOR THE NEW 
C ITERATION STEP. 

DC 18 1=1 ,M 

It WORK! I I = WORK! I I/R 
C 

C THE COMPUTATION OF the RESIOUALS AND COMPARISON OF THE 
C RF SI DUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE 
C ITERATION. IF THE INFINITE NORM OF THE RESIDUAL VECTOR IS 
C GREATER THAN THE INFINITE NORM OF THE PREVIOUS RESIDUAL 
C VECTOR THE COMPUTED EIGENVECTOR CF THE PREVIOUS STEP IS 
C TAKEN AS THE FINAL EIGENVECTOR. 
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R1 = 0.00 
CO 18 1=1, * 

T * 0.00 
00 17 J=t ,N 

17 t = T*A( ! , J I *W0RK( J I 

r * mm 

IFIRl.GE. TIGO TO 18 
Rl= T 

18 CONTINUE 

IF! ITER.EQ. 1>G0 TO 19. 

IFIPREVIS.LE.RI I GO TO 2A 

19 00 20 I =1 ,M 

20 VECRII.IVEC) = W0RKU 1 
FREVIS = R1 
IFINS.EQ.l )GO TO, 2<t 

IFI ITER .GT.6IG0 TO 25 
ITER = ITERU 
IFIR .LT.flOUNOIGO TO 21 
AS = l 

GAUSSIAN ELIMINATION OF THE RIGHT-HAND SIOE VECTOR 

21 K = M-l 

00 23 I =1 , K 
R = WORK! I ♦ 1 I 
IFI IHORKUI .EO.O ICO TO 22 

WORKlt + tl =WORK ( I I ♦NORM (I*l)*A( (♦1,11 
WORM I ) = R 
GC TO 23 

22 WORK 1 1 ♦ 1 ) *NORK! U1 J+WORKI l»*A( !♦! , It 
2? CONTINUF 

GC TC 12 

2A INC IC ( IVEC) = 2 

25 IF(M.EO.N)GO TO 27 
J = M* 1 

CO 26 I=J,N 

26 VECR ( I ,IVEC) = 0.00 

27 RETURN 
EKO 
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SUBROUTINE COMPVEIN.NM,M,l VEC ,A , VEC R ,H , E VR .E VI ,INDIC, 

1 I WORK, SUBDI A .WORM ,WCR K2 .WORK, E PS , E X I 
INTEGER I ,I1,I2,ITER.1VEC,J.K»L,M,N,NM,NS 
C IMFNS ION AINM, 1I.VECRI NM, 1),HCNM,1I ,EVR INM) .EVKNMI, 

1 INDICINMI , 1 WORK ( NMI , SUBOIAINM) t WORK 1 ( NM ) , WORK 2 ( NM ) , 
2WORKINM) 

C 

C THIS SUBROUTINE FINOS THE COMPLEX EIGENVECTOR OF THE REAL 
C UPPER-HESSENBERG MATRIX OF ORDER' N CORRESPONDING TO THF 
C COMPLEX EIGENVALUE WITH THE REAL PART IN EVRUVEO AND THE 
C CORRESPONDING IMAGINARY PART IN EVIIIVEC). THE INVERSE 
C ITERATION METHOD IS USED MODIFIED TO AVOID THE USE OF 
C CCMPLEX ARITHMETIC. 

C THE MATRIX ON WHICH THE INVERSE ITERATION IS PERFORMED IS 
C BUILT UP IN THE ARRAY A BY USING THE UPPER-HESSENBERG 
C MATRIX PRESERVED IN THE UPPER HALF CF THE ARRAY H AN C IN 
C THE ARRAY SUBOIA . 

C NM DEFINES THE FIRST DIMENSION OF THE TWO. DIMENSIONAL 
C ARRAYS A , VECR AND H. NM MUST BE EQUAL TO CR GREATER 
C THAN N. v 

C M IS THE CROER OF THF SUBMATRIX OPT A INEO BY \ A SUITABLE 
C DECOMPOSITION OF THE UPPER-HESSENBERG MATRIXsIF SOME 
C SUBD IACONAL ELEMENTS ARE EQUAL TO ZERO. THF VALUE Of " 15 
C CHCSEN SO THAT THE LAST N-M COMPONENTS OF THE COMPLEX 
C EIGENVECTOR ARE ZERO. 

C 

C THE REAL PARTS OF THE FIRST M COMPONENTS OF THE COMPUTED 
C COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST M PLACES OF 
C THE COLUMN WHOSE TOP ELEMENT IS VECRIl.IVEC) AND THE 
C CCRRFSPCNDING IMAGINARY PARTS OF THE FIRST M COMPONENTS OF 
C THE COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST « 

C PLACES CF THE COLUMN WHOSE TOP ELEMENT IS VECRI It IVFC- 1) . 

C 

C THF ARRAY I NO IC INDICATES THE SUCCESS OF THE ROUTINE AS 
C FCLLCWS 

C VALUE CE INOICIII EIGENVECTOR I 

C l NOT FOUND 

C 2 FOUND 

C THE ARRAYS I WOR K ,WORKl , WCRK2 ANC WORK ARE THE WORKING 
C STORES USED OUR ING THE INVERSE ITERATION PROCESS. 

C EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS 
C ZERO IN THE PROGRAM. EPS = IEUCLIOIAN NCRM OF H)*EX, WHERE 
C EX = ?**(-T ) . T IS THF NUMBER OF BINARY DIGITS IN THE 
C MAMTISSA OF A FLCAT I KG PCINT NUMBER. 

C 

EKS I = E VR ( IV EC I 
ETA = E VI ( IVECI 

C THE MODIFICATION OF THE EIGENVALUE IFKSI ♦ !*ETA» IF MORE 
C EIGENVALUES ARE EQUAL. 

IF ( I VEC.EQ. MJGC TC 2 
K = IVEC+l 
P = 0.00 
DO 1 I=K,M 

IFIFKSI .NE.EVRI I HGO TO l 
I F ( ABSIETAI.NE. ABSIEVMIJHGO TO I 
R = R ♦ 3.00 
1 CONTINUE 
P = R*EX 
EKSI = FK SI *R 
ETA = ETA »R 

THE MATRIX I (H-FKS!*M*IH-FK5I0! ) ♦ ! ET A*ET A )• I ) IS 
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C STOREO INTO THE ARRAY .A- 

2 R = FKS I*FKS ! ♦ ET A*ETA 
S = 2. OC*FK SI 
t = M-l 
CC 5 1 = 1, M 
OC 4 J=! ,H 
0 = 0.00 



A ( J , I ) 

= 0.00 

-■ 

DO 3 K 

= I,J 

3 

0 = 0 *H( IfK )*H(K,J ) 

4 

All .J) 

= 0-S*H( I , J) 

5 

All,!) = 

A(I ,t l«R 


CC 9 1=1, L 

R = SUBOI AIM 
A ( I ♦ l , I J = -S*R 

11 =1*1 

OC 6 J=1,I1 

6 A(J,I» = A< J,I )*R*H( J,I*1I 
IF ( I. EO. I ICO TO 7 

Aim, 1-1) = R*SUEC1 A ( 1-1 ) 

7 00 8 J=t,M 

8 A<T*l,J) = A ( I *1 » J I *R*H( 1,4 I 

9 CONTINUE 
C 

C THE GAUSSIAN ELIMINATION OF THE MATRIX 

C 1 1H-FKSI*! )*(H-FKSl*I I ♦ 1ETA*ETAI*I) IN THE ARRAY A. THE 
C ROW INTERCHANGES THAT OCCUR ARE INDICATED IN THE ARRAY 
C I WORK. ALL THE MULTI PLIERS ARE STORED IN THE FIRST AND IN 
C THE SECOND SUBO I AGONAL OF THE ARRAY A. 

K = P-1 
OC 18 I=1,K 
11=1*1 

12 = 1*2 
I WORK ( 1 1 = 0 
I F ( I .EQ.K I GO TO 10 
IF(A(I*2, II. NE. 0.001 GO TO 11 
IF(A(I*1,I).NE. 0. 00 1 GO TO 11 
IF1A1 I, I) .NE.O.COIGO TO 18 
A ( 1 , 1 ) = EPS 
GO TO 18 

IFU.EQ.KIGC TC 12 

I E ( ABS(AU*1 ,! M.GE. ABS ( All *2 , 1 III GO TC 12 
I F ( ABS< A< I, I I) .GE . ABS1AI 1*2,1 II I GO TO 16 
L = 1*2 
I WORK! 1 1 = 2 
GO TO 13 

IF( ABSIAII ,1) I.GE. ABS I A'( I.*l , I 1 1 I GO TO 15 
L = 1*1 
IWORKII) = 1 

00 14 J=I ,M 
R = A ( I , J ) 

AH , Jl = A( L , J) 

AIL , J I = R 
IF ( I .NE.K I GO TO 16 
12 = II 
DO 17 L= I 1* 12 

R = -AIL, D/AIT, I) 

AIL, I) = R 
DO 17 J=Il,M 

AU,J) = A( L , Jl *R* A 1 1 , J I 


1C 


C 

11 


12 


C 

13 

14 

15 

16 
17 
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18 CONTINUE 

IF(A|M,M).NE.O.OO)GO TO 19 
AH', Ml = EPS 
C 

C THE VECTOR (1,1,. ...11 IS STORED INTO THE RIGHT-HAND SIOE 
C VECTORS VECR ( ,IVEC1 AND VECRI .IVEC-ll REPRESENTING THE 
C CCMPLEX RIGHT -HANC SICE VECTOR, 

19 DO 21 1=1 ,N 

IFII.GT.MIGO TO 20 ' • 



VECRII, IVEC 1 = 

• 

o 

o 


VECRII, IVFC-ll 

= 1.00 


GO TO 2 1 


20 

VECRII, IVEC1 = 

0 .00 


VECRI I ,1 VEC-1 1 

= 0.00 

21 

CONTINUE 



C THE INVERSE ITERATION IS PEPFORMEC ON THE MATRIX UNTIL THE 
C INFINITE NORM OF THE RIGHT-HAND SIDE VECTOR IS GREATER 
C THAN THE 8CUN0 OEFINEC AS O.Ol/IN*EXl. 

BOUND = 0.0100/IEX* F LOAT (Nil 
NS = 0 
ITER =1 
00 22 I=l,M 

22 WCRKI I 1 = HI I, Il-FKSI 
C 

C THE SEQUENCE OF THE COMPLEX VECTORS ZISI = PISIM*3IS» AND 
C W (S -»1 1= UIS + 1 1H*V(S*1 1 IS GIVEN BV THE RELATIONS 
C (A - <FKSl-!*ETA)*I l*MS + l I = ZISI AND 
C Z f S-H 1 1 = WISMI/MAXI W( S*l)l. 

C THE FINAL WISI IS TAKEN AS THE COMPUTED EIGENVECTOR. 

C 

C THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR 
C «A-FKS!*I1*P(SI-ETA*Q(S1. A IS AN UPPER-HESSENBERG MATRIX. 
22 00 27 1=1, M 

0 = WORK ( I I HVECR ( I, IVEC 1 
IFII.EQ.il GO TO 24 
D = OhSURDIAI I-1I*VECR(I-1,IVECI 

24 L = 1*1 
IFIL.GT.MIGC TC 26 
DO 25 K=L,M 

25 0 * OHH( I,K 1*VECR1K, IVEC 1 

26 VECRII, IVEC-ll = O-ET A* VECR 1 1 , IVEC- 1 I 

27 CONTINUE 
C 

C GAUSSIAN ELIMINATION OF THE RIGHT-HAND SIDE VECTOR. 

K = M-l 
CC 28 1=1 »K 

L = I ♦ I WORK ( I 1 
R = VECRIL, IVEC-ll 
VECRIL ,IVEC-l 1 = VECRI I, IVEC-l I 
VECRI I , IVEC-ll = R 

VECR I I ♦ 1 , IVEC-l 1 = VEC R I I ♦ 1, I VEC- 1 1 *A I I»l,II*R 
IFII.EO.KIGC TC 2fl 

VECRI 1*2, IVEC-ll = VECRII *2,1 VEC-1 1*A( 1*2,11*8 

28 CONTINUE 

THE COMPUTATION OF THE REAL PART UIS*ll OF THE COMPLEX 
VECTOR WIS+ll. THE VECTOR U(S*l» IS O0TAINEO AFTER THE 
BACK SUB STITUTl ON. 

CO 31 l=l,M 
J = M-t*l 

0 = VECRI J, IVEC-ll 
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IMJ.EC.MIGO TO 30 
L = JM 
00 29 K=L , M 
01 = AIJ.KI 

29 0 = 0-01*VFCR1K, IVEC-1 I 

90 VECRU, IVEC-1) = D/AIJ.J) 

31 CONTINUE 
C 

C THE COMPUTATION OF THE IMAGINARY PART V(S*1) OF THE VECTOR 
C WISU >, WHERE VISU) = ( PI S I-U-FK S I *11 *UI S* 1 II /E TA . 

00 35 1=1, H 

0 = WORKU )*VECR< I , IVEC-1! 

IFII.EQ.l ICC TC 32 
0 = D+SU8DIAI 1 — 1 1 * VECR ( I -l » I VEC-l I 

32 L = 1*1 
IFIL.GT.MIGO TC 34 
no 33 K=L , M 

33 C = 0*H( I, K 1*VECR0C, IVEC-1) 

34 VECRU, IVEC) = IVECRt I , IVEO-OI/ETA 

35 CONTINUE 
(C 

C THE COMPUTATION CF UNFIN. NORM OF W(SM))**2 . 

L = l 
S = 0.00 
00 36 I =1 , M 

Wt = VECRU, IVEC) 

H? = VECRU. IVEC-1 ) 

R = W1*W1 ♦ W2*H2 
IFIR.LE.SIGO TO 36 
S = R 
L = I 

36 CONTINUE 

C THE COMPUTATION OF THE VECTCR 2<S*l), WHERE MS*1>* W(SM)/ 
C (COMPONENT OF W(S*1) WITH THE LARGEST ABSOLUTE VALUE) . 

U = VECRIL, IVEC-1) 

V = VECR ( L , IVEC I 
00 37 1=1, M 

B = VECRU, IVEC) 

R = VECRU .IVEC-1 ) 

VECRU, IVEC) = <R*U ♦ B*V)/S 

37 VECRU , IVEC-1 ) = <B*U-R*V)/S 

C THE COMPUTATION OF THE RESIOUALS AND COMPARISON CF THE 
C RESICUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE 
C ITERATION. IF THE INFINITE NORM OF THE RESIDUAL VECTOR IS 
C GREATER THAN THF INFINITE NORM OF THE PREVIOUS RESIDUAL 
C VECTOR THE COMPUTED VECTOR OF THE PREVIOUS STEP IS TAKEN 
C AS THE COMPUTED APPROXIMATION TO THE EIGENVECTOR. 

B = 0.00 

CC 41 1=1, M 

R = WORKU)*VECRU , IVEC-1) - ET A* VEC R (I , IVEC ) 

U = WORK I I )*VECR ( I , IVEC ) ♦ ETA*VECR U , IVEC-1 ) 
IFII.EQ.l ) GO TC 38 
R = R»SUBD!A( I -l ) * VECR (I - l , IVEC-1) 

U = U+SUBDI A( 1-1 )*VECR ( I- 1, IVEC ) 

3 6 1 * IM 

IF(L.GT.M)GO TO 40 
OC 39 J=L,M 

R = R*H(I ,J)*VECR< J.IVEC-1) 

39 U = U*H( I,J)*VECR(J. I VFC ) 

40 U = R*R ♦ U*U 
IFIB.GE.U)GO TO 41 
B = U 
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*1 CONTINUE 

I FI ITER .EQ.l )GO TO 42 
IFIPREVIS.LE.BIGC TO 44 
*2 00 43 1=1 ,N 

WORK 1 1 !> = VECRI I t 1VEC 1 

43 WCRK2I!) = VECR I 1 , IVEC-l ) 

PREVIS = B \ 

tFINS.EQ'.llGO TO 46'- 

IFI ITEP «GT-.6 IGC TC 47 s " 

ITER = I TER* 1 

IF( BOUND. GT . SQRTISIIGO TO 23 
KS = l 
GO TO 23 

44 00 45 I =1 .N 

VECRI I , IVECI = WORK 1 1 I > 

45 VECRI I. IVEC-i »=W0RK2I 1 1 
44 I NO I C I f VEC-1J = 2 

INCICt IVECI = 2 
47 RETURN 
END 
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SUBROUTINE RRAY 


C 

C 

C 

C 

C 

c 

c 

c 

f. 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

r. 

r 

c 

c 

c 

c 

r 

C 

C 

C 

C 

C 

C 

C 

C 

r. 

f. 

c 

r 

c 

c 

r. 

c 

c 

c 

e 

c 

c 

r 


PURPOSE 

CONVERT DATA ARRAY PROP SINGLE TO COUBLE DIMENSION OR VICE 
VERSA. THIS SUBROUTINE IS USED TO LINK THE USER PROGRAM 
WHICH HAS DOUBLE DIMENSION ARRAYS AND THE SSP SUBROUTINES 
WHICH OPERATE ON ARRAYS OF OA TA IN A VECTOR FASHION. 

USAGE 

CALL ARRAY < MODE , I , J, N , M, S ,D> 

DESCRIPTION OF PARAMETERS 

MODE - CODE INDICATING TYPE OF CONVERSION 

1 - FROM SINGLE TO OOUBt E DIMENSION 
?. - FROM DOUBLE TO SINGLE DIMENSION 
I - NUMBER OF ROWS IN ACTUAL DATA MATRIX 

.1 - NUMBER CF COLUMNS IN ACTUAL DATA MATRIX 

N - NUMBFR OF ROWS SPECIFIED FOR THE MATRIX D IN 

DIMENSION STATEMENT 

M - NUMBER CF COLUMNS SPECIFIED FOR THE MATRIX D IN 
DIMENSION STATEMENT 

S - IF MODE* 1 , THIS VECTOR IS INPUT WHICH CONTAINS THE 

ELEMENTS OF A OATA MATRIX OF SIZE I BY J. COLUMN 1*1 
Of DATA MATRIX FOLLOWS COLUMN I, ETC. IF MOOE=2 , 

THIS VECTOR IS OUTPUT REPRESENTING A DATA MATRIX OF 
SIZE I BY J CONTAINING ITS CCLUMNS CONSECUTIVELY. 

THE LENGTH OF S IS TJ. WHERE IJ=I*J. 

0 - IF MOOF=l, THIS MATRIX OF SIZE N BY M IS OUTPUT, 

CONTAINING A OATA MATRIX OF SIZE I BY J IN THE FIRST 
I ROWS ANO J COLUMNS. IF MOOE=2, THIS N BY M MATRIX 
IS INPUT CONTAINING A DATA MATRIX OF SIZE I BY J IN 
THE FIRST I ROWS AND J COLUMNS. 

REMARKS 

VECTOR S CAN BE IN THE SAME LOCATION AS MATRIX D. VECTOR S 
IS REFFRREC AS A MATRIX IN OTHER SSP ROUT IMES , . S INCE IT 
CONTAINS A OATA MATRIX. 

THIS SUBROUTINE CONVERTS ONLY GENERAL DATA MATRICES (STORAGE 
MCOE OF OI. 

SUBROUTINES ANO FUNCTION SUBROUTINES REQUIRED 
NONE 

METHOD 

REFER TO THE OtSCUSSION ON VARIABLE DATA SIZE IN THE SECTION 
DESCRIBING OVERALL RULES FOR USAGE IN THIS “ANUAL . 


C 

r 

C 

r 

r 

c 

r 


SUBROUTINE RRAY I MOOE, I , J, N, M, S ,0 > 

01 ME NS! ON Sill ,0(11 

N I* A - I 

TEST TYPE OF CONVERSION 
IF I MC.OF - II 100, 100, 120 

CONVERT FROM SINGLE TO OOUBLE DIMENSION 
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ICO IJ=I*J* 1 

oo no k*i • j 

NM=NM-NI 

cc no l»i ,1 

IJ*IJ-1 

NM=NM-1 

no c i n** > =s ( i j) 

GO TO 140 
C 

C CCNVERT FROM COUBLE TO SINGLE DIMENSION 

C 

1 2C IJ=0 
M <=0 

no no k=i,j 

00 125 1=1, I 
IJ* 1J*1 
M* = M**1 

125 SI IJ 1=01 NM 1 
130 M»=NN*N! 

C 

140 RETURN 
ENO 
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SUBROUTINE MPRC 
PURPOSE 

MULTIPLY TWO GENERAL MATRICES TO FORM A RESULTANT GENERAL 
MATRIX 

USAGE 

CALL MPRDI A,B ,A,N,M,L! 

DESCRIPTION OF PARAMETERS 

A - NAME OF FIRST INPUT MATRIX 
8 - NAME OF SECOND INPUT MATRIX 
R - NAME OF OUTPUT MATRIX 
N - NUMBER OF ROWS IN A 

M - NUMBER OF COLUMNS IN A AND ROWS IN 6 

L - NUMBER CF COLUMNS IN B 

REMARKS 

ALL MATRICES MUST BE STOREO AS GENERAL MATRICES 
MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX A 

MATRIX R CANNOT BE IN THE SAME LOCATION AS MATRIX B 

NUMBER OE COLUMNS OF MATRIX A MUST BE EQUAL TO NUMBER OF ROW 
OF MATR IX B 

SUBROUTINES ANO FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

THE M BY L MATRIX B IS PRE MUL TI PI I ED ft Y THE N BY M MATRIX A 
ANC THE RESULT IS STOREO IN THE N RY |. MATRIX K. 


SUBROUTINE MPRO I A ,B ,R ,N,M, LI 
CIMENSION A(1),B(1),R(1| 

I R = C 
!K=-M 

CC 10 K=1,L 
I K= I K+M 
CO 10 J= 1,N 
IR=IR*1 
J I = J-N 
IB= IK 
P 1 1 R I =0 
00 1C I=1,M 
Jl=J 1 *N 
I B= I B*1 

1C Rl IRI=RUR)*A(JII*B( IB) 

RETURN 

END 
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SLBROLTl NE I NV 


PURPCSE 

INVERT A MATRl X 


USAGE 

CALL INV <A,N,L,M) 

DESCRIPTION OF PARAMETERS 

A - INPUT MATRIX, DESTROYED IN COMPUTATION AND REPLACED BY 
RESULTANT INVERSE. 

N - OROER CF MATRIX A 
D - RESULTANT DETERMINANT 
L - WORK VECTOR OF LENGTH N 
m - WORK VECTOR OF LENGTH N 

REMARKS 

MATRIX A MUST BE A GENERAL MATRIX 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOO 

THE STANOARO GAUSS-JORDAN METHOD IS USED. THE DETERMINANT 
IS ALSO CALCULATED. A DETERMINANT OF ZERO INDICATES t H AT 
THE MATRIX IS SINGULAR. 


SUB ROUT INF INV (A.N.L.MJ 
DIMENSION All), Lilt, Mill 

C 

C 

C 

C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE 

C C IN COLUMN 1 SHOULD BE REMOVED FROM THE CP'JBLE PRECISION 

C STATEMENT WHICH FOLLOWS. 

C 

C DOUBLE PRECISION A ,D ,BI GA .HOLD 

C 

C THE C MUST ALSO BE REMOVEO FROM DOUBLE PRECISION STATEMENTS 

C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS 

C ROUTINE. 

r 

C. THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO 

C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. ABS IN STATEMENT 

C 10 MUST BE CHANGEC TO DABS. 

C 

C 

c 

C SEARCH FOR LARGEST ELEMENT 

C 

0 = 1.0 

NK=-N 

CO 80 K=l,N 
NK=NK*N 
L(K)=K 
MIK )=K 
KK=NK*K 
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B! GA =A( KK) 
rn 20 J*K,N 

iz=n*ij-i> 

00 20 I=K,N 

1 J= I 2 ♦ I 

1C IF( ARSI BIGA) - ABS( A( IJI ) > 15,20,20 
15 eiGA=A( IJ) 
l IK 1 * f 
I K J =J 

2C CONTINUE 

INTERCHANGE RONS 
J = L <K I 

IFIJ-KI 35,35,25 
25 KI=K-N 

CO 30 1=1, N 
K I =K l * N 
HOLP =- A I K I • 

J I = K I-KfJ 
A I K l ) = A(J1» 

2C A ( J I ) = H0L0 

INTERCHANGE columns 

35 I = M ( K I 

IF(I-K) *5,45, 3B 
2f JP=N* ( I -1 ) 

CO AO J=l,N 
JK=NK*J 
JI=JP*J 
FOL C=- A I JK ) 

AUK) = MJI) 

AC A ( J I ) =MC!L0 

C1VIDE COLUMN BY MINUS PIVOT (VALUE OF <>IV0T ELEMENT IS 
CONTAINED IN B1GA1 

A 5 IFIBIGAI AB.A6.A8 
Af D=C.C 
RETURN 

AB CC 55 1=1, N 

IF(I-K) 50,55,50 
50 IK= NK ♦ I 

A( IK MAI IK1/ (-RIGA) 

5= CONTINUE 

RFOUCE MATRIX 

00 65 1= 1, N 
IK=NK*I 
HQLO=A ( IK) 

1 J= 1-N 

CC 65 J=l,N 
I J=I JFN 

I F ( I-K I 60, 65, 60 
60 IF(J-K) 62,65,62 
62 KJ=IJ-I*K 

A ( IJ J=K)tO*A(KJ MAI IJ) 

65 CONTINUE 

OIVIDE ROW BY PIVOT 
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KJ=K-N 
CC 75 J=1,N 
K J = K J«-K 

IFCJ-K1 70,75,70 
70 ACKJ)=ACKJ)/e!CA 
75 CONTINUE 

PROOUCT OF PIVOTS 

0=0*8 1 GA 

REPLACE PIVOT ev R EC tPROCAL 

A ( KK >= 1.0/B IGA 
80 CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 


K=N 

ICC K = (K-ll 

IFIKI 150,150,105 
1C5 I =1 ( Kl 

IF! I-K I 120,120, 1 C 8 
108 JC=N*(K-ll 
JR = N* ( I -1) 

ro no j=i,n 

JK=JC*J 
HOLD = A ( JK1 
J I = JR ♦ J 
AC JK ) = - A ( J I I 
110 ACJII =HOL 0 
12C J=M C K I 

IHJ-KI 100,100,125 
125 Kl=K-N 

no 13" 1=1, N 
K 1 = K I*N 
HC(.0 = A( K [ I 
Jl=K!-K«-J 
A ( K l ) = - A ( J I | 

13C ACJII =HOL D 
GO TO 100 
150 RETURN 
END 


79 



FUNCTION RSHL(T) 

A = 14.39 
B = 9.35 

EPS=1.0-(B/A)*(B/A) 

RSHL=B*B/( A*< 1.0-EPS* COS(T)* COS ( T ) ) **1 . 5 ) 

RETURN 

END 


FUNCTION RRRT ( T) 

A = 14.39 
B = 9.35 

EPS=1.0~(B/A)*(B/A) 

RRRT= SORT < 1 .O-EPS* COS ( T ) * COS ( T ) ) *A*EPS* S IN < 2 • 0*T ) *1 . 5/ ( B#B ) 

RETURN 

END 


FUNCTION RSHLT(T) 

A = 14.39 
B = 9.35 

EPS=1.0-(B/A>* (B/A> 

RSHLT=-1. 5*B*B*EPS* S J N < 2 . 0*T ) / ! A* ( 1 .O-EPS* COS(T)* COS ( T ) ) *■*■■?. » 5 ) 

RETURN 

END 


) 
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APPENDIX C 


DICTIONARY OF VARIABLES USED IN THE MAIN PROGRAM 


AA 

Length of the shell 

ABN 

nn 

ABNA 

n 3 n 

ABNB 

nn a 

ABN2 

a - a 
n n 

AK(MN3,MN3) 

Stiffness matrix, dimension should be = KRRR 

AM(MN3,MN3) 

Mass matrix, dimension should be = KRRR 

AN 

n 

AN2 

n a 

AR(K) 

th 

Cross-sectional area of the k kind or ring, dimension 
should be i E 

AS(L) 

Cross-sectional area of the kind of stringer, 

dimension should be s NL 

BC(2,4) 

Data block defining various boundary conditions, 
dimension should be 2 * 4 

BCR(2) 

The name of the boundary condition read-in, dimension 
should be 2 

BN 

n 

BN2 

i; 3 

C(8) 

Temporary work vector used in the stringer equations, 
dimension should be 8 
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1 


cc 

CG 

CN 

CNB 

CR(K,40) 

CS 

D 

DR(9) 

DRV ( 5 ) 
E1R(K) 

E2R(K) 

ElRK 

E2RK 

EC 

ER(K) 

ES(L) 

EVI(MN3) 


CN X CNB 


Temporary work variable 

Cos (TN) for NSA = 0; Sin (TN) for NSA = 

Cos (TNB) for NSA = 0; Sin (THB) for NSA = 1 

C^ to constants used in the ring equations (C3), 

and are defined in Appendix D of Reference 1, The first 
dimension should be ^ NK, and the second should be 40. 

CN X SNB 


Isotropic plate flexural stiffness 
Vector of the integrands of the circumferential 
integrals ISl^ to ISl^, dimension should be 9 
Vector of the integrands of the circumferential 
integrals 182.^ to IS2,. , dimension should be 5 
z-distance of the shear center of the k fc ^ kind of ring 
from the middle surface of the shell, dimension should 
be NK 

z-distance of the centroid of the kind of ring from 
its shear center, dimension should be s NK 


'Irk 


2rk 

Shell Young's Modulus 

Young's modulus of k^ kind of ring. 


dimension should 


be ^ NK 

Young's modulus of kind of stringer, dimension 
should be s NL 

Vector of imaginary part of the eigenvalues, dimension 
should be = KRRR 
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EVR(MN3) 

GJR(K) 

GJS(L) 

H 

I 

IBC 

IM 

IN 

INDIC (MN3) 


INN 

IR 

ITEMP,IY1) 

IY2, IZ1, IZ2) 
J 

JBC 

JN 

JNN 


Vector of real part of the eigenvalues, dimension 
should be = KRRR 

The torsional stiffness of the kind of ring, 
dimension should be > NK 

The torsional stiffness of the kind of stringer, 
dimension should be £ NL 
Shell thickness 

Row index of [a] , [d], [e], [n], [NN], and [p] sub- 
matrices 

Temporary work variable 
Temporary work variable 

Row index of [b], [f], [q], and [r] submatrices 
This array indicates the success of the subroutine 


EIGENP as 

follows : 


INDIC(I) 

EIGENVALUE I 

EIGNEVECTOR I 

0 

not found 

not found 

1 

found 

not found 

2 

found 

found 


Row index of [c] and [ s] submatrices 
(see the listing of the program) 

Temporary work variables 

Column index of [a] , and [n] submatrices 
Temporary work variable 

Column index of [dJ , [b], [nn] , and [q] submatrices 
Column index of [e], [f], [c] , [p], [r] , and [s] 
submatrices 
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JTEMP 

KG.KK 

KQ 

KRRR 

LC(MN3) 

LL 

MC (MN3) 

MD 

MMAX, MMIN 

MN3 

MS 

MSA 

NBC 


NCHNG, ND 
NDC.NEC 
NEO 
NEIXT 


Temporary work variable 

(see the listing of the program) 

Temporary work variable 

Dimension of [AK] , [AM] , [VECR] , { EVR} , { EVl} , { INDlc} , 
{ Lc} , and {mc} > MN3 

Temporary vector used by INV (matrix inversion) sub- 
routine, dimension should be = KRRR 
(See listing of the program) 

Temporary vector used by TNV (matric inversion) sub- 
routine, dimension should be = KRRR 
Temporary work variable 
(See the listing of the program) 

Order of the mass, stiffness, and modal matrices 

Total number of axial mode components considered in the 

displacement series 

(See the listing of the program) 

The code number assigned for different boundary 
conditions as follows: 

1 for clamped-free 

2 for freely supported 

3 for clamped-c lamped 

4 for free-free 

Temporary work variables 

(See the listing of the program) 

Temporary work variable 
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NG, NK, NL 
NMAX, NMIN 
NNK (K) 

NNL(L) 

NNR(I) 

NQUIT 


NR(NK ,NK) 
NS 

NSA.NWEV, 

NWK , NWM 

PC 

PHI 

PI 

PI2 

PR(K) 

PS (L) 

R(9) 


(See the listing of the program) 

Vector of the number of rings of kind of ring, 
dimension should be > NK 

Vector of the number of stringers of kind of 
stringer, dimension should be > NL 

Temporary vector containing centroidal information of 
different kinds of rings, dimension should be ^ NK 
1 in the 80th column of a blank card, which when placed 
at the end of sets of data, signifies the end of the 
data 

Temporary vector containing centroidal information of 
the rings, dimension should be NK * NK 
Total number of circumferential mode components con- 
sidered in the displacement series 

(See the listing of the program) 

Mass density of shell 
Temporary work variable 
it = 3.14159 

2tt 

Mass density of k*"* 1 kind of ring, dimension should 
be 5 NK 

th 

Mass density of H kind of stringer, dimension 
should be 2 NL 

Vector of the circumferential integrals ISl^ to ISl^, 
dimension should be = 9 
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Rl(8) , R2( 10) , 
R3(2), R4(5) , 
R5( 18) , R6(ll) 


RCG(K) 


RI(K ,54) 


RING1 to RING6 

RR1(8), RR2(10) 
RR3(2), RR4 (5) 
RR5(18), and 
RR6(11) 

RRRT (0) 

RSHL(0 ) 

RV(5) 

RX(K, I) 


The vectors, of the integrands of the circumferential 
integrals of the ring, IR^-IRlg, , IR2 1 -IR2 1() , IR3 1 
and IR3,,, IR4 -IR4_, IR5 -IR5 10 , and IR6 -IR6.,, 
respectively; dimensions should be,= 8, 10, 2, 5, 

18, and. 11, respectively 

Vector of centroidal distances of various kinds of 
rings, dimension should be > NK 

Temporary work vector for saving the 54 ring integrals 
IR1 to IR6^^. The first dimension should be 2 NK, and 
the second dimension should be = 54 

Subroutines defining the integrands of the circumferen- 
tial integrals of the ring, IRl^ to IR6^ 

The vectors of the circumferential integrals of the 
ring, IR1, - IR1 Q , IR2, - IR2 in , IR3, and IR3„, IR4, - 

1 o 1 10 1 2 . 1 

IR4 , IR5 - IR5 Q , and IR6 - IR6 , respectively; 

-> 1 ,^lo * H 

dimensions should be equal to 8, 10, 2, 5, 18 and 11, 
respectively. 

Function subroutine furnished by the user of the 

program to evaluate (1/R),g at a given value of 0 

Function subroutine furnished by the user of the 

program to evaluate R at a given value of 6 

Vector of circumferential integrals of the shell, 

IS2^ to IS2,., dimension should be 5 

th 

Array of the x-locations of the k kind of rings, 
the first dimension should be 2 NK, and the second 
dimensions should be 2 the largest element of the vector 


NNK(K) 



SI to S8 
SC 

SHELL 1, SHELL 2 

SN 

SNB 

SR 

SR2 

SS( 1 ,30) 

SSS 

ST(75) 

SUM( 18) 

T(L,I) 


TITLE1(7) 

TITLE2(7) 

TN 

TNB 


S. to S 0 , constants used in Eqs. (Cl), and are defined 
1 o 

in Appendix D of Volume I 
SN X CNB 

Subroutines defining the integrands of the circumferen- 
tial integrals ISl^ to IS2^ of the shell equations 
Sin (TN) for NSA = 0; Cos (TN) for NSA = 1 
Sin (TNB) for NSA = 0; Cos (TNB) for NSA = 1 
Radius of the shell, R 
R 3 

SS^ to SS 3Q , constants used in Eqs. (C2), and are 
defined in Appendix D of Volume I, The first dimension 
should be 2 NL, and the second dimension should be = 30 
SN X SNB 

Intermediate terms of the stringer equations, dimension 
should be = 75 

Temporary work vector used by GAUSS (numerical inte- 
gration) subroutine, dimension should be = 18 
Array of the 0- locations of the SL ^ kind of stringers, 
the first dimension should be 2 NL, and the second 
dimension should be £ the largest element of the vector 
NNL(L) 

Title of the run 
Title of the run continued 
n X t (L,I) 

— ii X T (L.I) 

\ 

\ 


I 
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TS(L,42) 


U 

VECR(MN3, MN3) 
X(5 ,IM) 


XI to X5 
XIR(K) 

XK 

XNU 

XR(1), XR(2) 

XXI, XX2 
XXX(2 , K, IM) 


XXXX(MN3 2 ) 


to constants used in Eqs. (C2) and are defined 

in Appendix D of Volume I, the first dimension should 
be £ NL, and the second dimension should be = 42 
Equal to number of binary digits in the mantissa of a 
double precision, floating point number 
Eigenvector (modal) matrix, dimension should be 
= KRRR X KRRR 

A temporary work matrix, containing the longitudinal 

integrals IX^ to IX,. for every combination of m and m; 

the first dimension should be = 5, and the second 

dimension should be = MS X (MS + l)/2 

IX^ to IX,. longitudinal integrals 

th 

The moment of the k kind of ring cross-sectional 
area about an axis parallel to X-axis passing through 
its centroid 

til 

X-location of the k ring 
Poisson's ratio 

Temporary work vector used for transferring X^ and 

values from XX subroutine to the main program 

X^, X^ (see eqs. (C8) in Appendix C of Volume I) 

A temporary storage three dimensional matrix 

containing the quantities X^ and X^ for every 

combination of m and m; the first dimension should be 

= 2, the second be ^ NK, and the third should be 

(MS -I- 1) MS 
2 

2 

Temporary work vector, dimension should be = (MN3 ) 



Y(MN3 2 ) 

Y1S(L) 


Y2S(L) 

YIS(L) 


YZIS(L) 

Z(MN3 2 ) 

Z1S(L) 

Z2S (L) 

ZERO 

ZIR(K) 

ZIS(L) 


2 

Temporary work vector, dimension should be = (MN3 ) 

y-distance of the shear center of the kind of 

stringer from the z-axis passing through its point 

of attachment, dimension should be ^ NL 

y-distance of the centroid of the kind of stringer 

from the shear center, dimension should be s NL 

The moment of inertia of the kind of 

stringer cross-sectional area about an axis parallel 

to y-axis passing through its centroid, dimension 

should be > NL 

th 

Product of inertia of the l kind of stringer cross- 
sectional area about y and z axes passing through its 
centroid, dimension should be s NL 

2 

Temporary work vector, dimension should be = (MN3) 

fch 

z-distance of the shear center of the & kind of 

stringer from the middle surface of the shell 

th 

z-distance of the centroid of the l kind of stringer 
from its shear center 

0.0, lower limit of the circumferential integrals < r 
shell and ring 

The moment of inertia of the k^ kind of ring cross- 
sectional area about z or z ' axes, dimension should 
be £ NK 

The moment of inertia of the kind of stringer 
cross-sectional area about an axis parallel to z-axis 
passing through its centroid, dimension should be ^ NL 
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APPENDIX D 

PREPARATION OF DATA FOR THE PROGRAM OF 
THE FREE VIBRATIONS OF RING- AND/OR 
STRINGER STIFFENED NONCIRCULAR 
CYLINDERS WITH ARBITRARY 
END CONDITIONS 




DATA 

No. of 
CARDS 

FORMAT 

ITEMS ON DATA CARD 

fj 

< 

(a) 

Name of the 

boundary 

condition 

1 

2A10, 

59x, 

11 

BCR, NQUIT 

pi 

w 

S3 

Id 

(b) 

General input 
parameters 

1 

2014 

NG, KG, LL, NL, KK, NK, NMIN, 
MMAX, MSA, NMIN, NMAX, MSA, 
NEW, IR, NWK, NWM, NWEV 

e> 

(c) 

Title of the 
run 

2 

7A10/ 

7A10 

TITLE 1, TITLE 2 

SHELL 

(a) 

Geometric and 
material pro- 
perties of 
shell 

1 

5E15.8 

PC, EC, XNU , H, AA 

Pi 

w 

o 

S3 

(a) 

Number of A ^ 
kind of 
stringers 

1 

14 

NNL(L) 

& 
c n 

(b) 

List of 8- 
locations of 
A fc b kind of 
stringers 

NNL(L) 

5 

5E15.8 

(T(L, I) , I - 1, NNL(L) ) 
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STRINGER 


DATA 


No. of 
CARDS 


FORMAT 


ITEMS ON DATA CARD 


(c) Geometric and 3 5E15.8 PS(L) , ES(L), AS(L) , 

material pro- Z1S(L), Z2S(L), Y1S(L), 

per ties of i th Y2S(L), ZIS(L), YIS(L), 

kind of stringer YZIS(L), GJS(L) 



(a) 

Number of k*"* 1 

1 

14 

NNK(K) 



kind of rings 





(b) 

List of X-lo- 

NNK(K) 

5E15.8 

(RX(K,I) , 1=1, NNK(K) ) 

o 


cations of k*-* 1 

5 



23 

M 


kind of rings 




04 

(c) 

Geometric and 

2 

5E15.8 

PR(K) , ER(K) , AR(K) , EIR(K) , 



material pro- 



E2R(K) , ZIR(K) , XIR(K) , 



perties of k fc b 



GJR(K) 



kind of ring 





Note: (1) All numerical data must be right justified. 


(2) No blank spaces are left before and in between data fields. 
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APPENDIX E 


COMPUTER OUTPUT 
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*** ******* ********************************************************* 


FREE VlbRATIONAL ANALYSIS OF STIFFENED GR UNSTIFFENED CIRCULAR OR 
NONCIRCULAR CYLINDERS MlTH ARBITRARY END CONDITIONS 

* ******** **** ****** ************* *********************************** 


GENERAL INPUT INFORMATION 


NG 

8 

KG 

4 

LL 

16 

NL = 

1 

KK 

11 

NK 

1 

MM I N = 

1 

MM AX = 

9 

MSA = 

1 

NMIN = 

1 

NMAX = 

1L 

NSA = 

0 

NEO = 

1 

IR = 

1 

NmK = 

0 

NWM = 

0 

Nw EV = 

0 








SHELL DATA 


SE WALL* S lb STRINGER AND 11 RING STIFFENED ELLIPTICAL 
CYL1NOER WITH A = 14.39 B = 9.35 


MASS OENSI TY 

MODULUS OF ELASTICITY 

POISSON* S RATIO 

THICKNESS 

LENGTH 

END CONDI TIONS 


0.2 5d80000U-03 LB :> EC. **2 / IN . **4 
0 • 10000000U 08 LB/IN.**2 
0. 3Q000000D 00 
0 ..>20000000-01 INCHES 
0.24000000D 02 INCHES 
FREELY SUPPORTED 
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STRINGER DAT A 


(THE UNITS ARE SAME AS THOSE OF SHELL OAT A ) 

TOTAL NUMBER OF STRINGERS = 16 

NUMBER OF DIFFERENT KINDS OF STRINGERS = I 


16 STRINGERS WITH THE FOLLOWING PROPERTIES 

MASS DENSITY = 0.25880C00D-03 MUD. OF ELAS. = 0. 10600000D 08 

AREA = 0.103687140 00 SHEAR CTR. (Zl)= -0.475000000-01 

SHEAR CTR. (Yl) = 0.0 CENTROID (Z2) = -0.233989590 00 

CENTROID (Y2J = 0.0 INERTIA I IZZ ) = 0 . 12 85071 70- 02 

INERTIA (IYYI = 0. 595710420-02 PROD. INER. ( IYZ)= 0.0 

TORSIONAL STIFFNESS = 0.91250000U 03 

LOCATED AT FOLLOWING THETA VALUES ( DEGREES J 

0.0 0.135000000 02 0.25500000D 02 0.518000000 02 

C.90CCOOOOD 02 0. 128200000 03 0.1^4500000 03 0.16630000D 03 

0.180000000 03 0.193500000 03 0.20550000D 03 0.23180000D 03 

0.270000000 03 0.30820000U 03 0.334500000 03 0.34650000U 03 


93 


KINO DATA 


(The units are same as those of shell data) 

TOTAL NUMBER OF RINGS = II 
NUMBER OF DIFFERENT KINDS OF RINGS = . 1 


11 RINGS WITH THE FOLLOWING PROPERTIES 

MASS DENSITY = 0.258800000-03 MOO. OF ElAST. = 0.10600000D 08 

AREA = 0. 10 36 U 7 140 00 SHEAR CTR. (El) = -0. <*75000000-0 1 

CENTROID (E2> =-0.233989590 00 INERTIA I 1ZZ > = 0. 128 507170-02 

INERTIA (IXX) = 0.5957 10420- 02 TORS. STIF.(GJ) = 0.91250000003 

LOCATED AT FOLLOWING X VALUES (INCHES) 

0.20000000D 01 0 .400000000 01 0.60000000D 01 0.80000000D 01 

0.100000000 02 0.120 00 00 OD 02 0.1<*OOOOOOD 02 0.180000000 02 

O.ldOOOOGOD 02 0. 200000000 02 0.22000000D 02 


94 


£1 GENV ALU bS IN HERTZ 


C. 52 5201 990 05 
0.443808860 05 
0 .325253081) 05 
C. 307913750 05 
0.270630280 05 
0.227547250 05 
0.197195290 05 
0.185293530 05 
0.1 7 7051850 05 
0.157811600 05 
0.144462940 05 
0.132438860 C5 
0.120629020 05 
0.103104730 05 
0.811092180 04 
0 .758435280 04 
0.6C2823540 04 
0.471157240 04 
0.740998740 03 
0. 17 02 8734D 04 
0.241361850 04 
0.284138540 04 
0.295882040 04 


0.460016450 

05 

0.38 7653 730 

05 

0.309578200 

05 

0. 291781 840 

05 

0.259449 520 

05 

0.233702 04 D 

05 

0. 198710190 

05 

0 .183604590 

05 

0.163912010 

05 

0.152548110 

05 

0 .145446830 

05 

0.136364 940 

05 

0.110746050 

05 

0.991873020 

04 

0. 792999230 

04 

0.618868420 

C4 

0.554182620 

04 

0.447110730 

04 

0.974276260 

03 

0.173906970 

04 

0.361588130 

04 

0.285624590 

04 

0. 310827990 

04 


0.486471890 05 
0. 332841530 C5 
0.313550720 05 
0.298557000 05 
0.250799860 05 
0.220058790 05 
0.192803320 05 
0. 178856710 05 
0.167839380 95 
0. 151812640 05 
0. 1427500 70 05 
0.127304990 05 
0. 107503390 05 
0.970742020 04 
0.781446350 04 
0.659350450 04 
0.497429840 04 
0.4075^4170 04 
0.115486790 04 
0 . 20087592 D 04 
0. 347187140 04 
0.292366990 04 


0.436157540 05 
0.340883090 05 
0.311164360 05 
0.267002370 05 
0.239849450 05 
0.22316/410 05 
0. 188148 71 D 05 
0. 172501550 05 
0.154738580 05 
0.147557350 05 
0.141590220 05 
0.13097o59D 05 
0.108287690 05 
0.986895570 04 
0.749336500 04 
0.600867310 04 
0.515371260 04 
0.437945790 03 
0. 133981750 04 
0 *37004882 U 04 
0.331427110 04 
0.311332820 04 
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